Lingfeng Duan

Chenglong Huang

Guoxing Chen

Lizhong Xiong

Qian Liu

Wanneng Yang

High-throughput estimation of yield for individual rice plant using multi-angle RGB imaging

Keywords: yield estimation, multi-angle imaging, individual rice plant, plant phenotyping, high-throughput

Modern breeding technologies are capable of producing hundreds of new varieties daily, so fast, simple and effective methods for screening valuable candidate plant materials are urgently needed. Final yield is a significant agricultural trait in rice breeding. In the screening and evaluation of the rice varieties, measuring and evaluating rice yield is essential. Conventional means of measuring rice yield mainly depend on manual determination, which is tedious, labor-intensive, subjective and error-prone, especially when large-scale plants were to be investigated. This paper presented an in vivo, automatic and high-throughput method to estimate the yield of individual pot-grown rice plant using multi-angle RGB imaging and image analysis. In this work, we demonstrated a new idea of estimating rice yield from projected panicle area, projected area of leaf and stem and fractal dimension. 5-fold cross validation showed that the predictive error was 7.45%. The constructed model achieved promising results on rice plants grown both in-door and out-door. The presented work has the potential of accelerating yield estimation and would be a promising impetus for plant phenomics.

Introduction

The growing population and the impact of changing climate on agriculture have put major crisis on world's food supply [START_REF] Rt | Plant phenomics: from gene to form and function[END_REF]. 70% more food is needed by 2050 according to the recent Declaration of World Summit on Food Security (www.fao.org/wsfs/world-summit/en/). As the staple food for over half of the world's population [START_REF] Yong | Deletion in a quantitative trait gene qPE9-1 associated with panicle erectness improves plant architecture during rice domestication[END_REF], rice production is significant for the food security of the world.

Modern breeding technologies are capable of producing hundreds of new varieties daily, so fast, simple and effective methods for screening valuable candidate plant materials are urgently needed [START_REF] Wanneng | High-throughput measurement of rice tillers using a conveyor equipped with x-ray computed tomography[END_REF].

Phenotyping the interested populations is widely deemed as the most laborious, costly and technically challenging part in plant breeding. Plant phenotyping has become one of the new bottlenecks in plant science and plant breeding [START_REF] Llorenç | Highthroughput phenotyping and genomic selection: The frontiers of crop breeding converge [J][END_REF][START_REF] Robert | Phenomics-technologies to relieve the phenotyping bottleneck[END_REF]. Plant phenomics enables the acceleration of progress in linking gene function, phenotype, and environmental responses, which will eventually accelerate plant breeding to develop new germplasm to meet the future food demanding.

Digital image analysis has been proved an effective tool in plant phenomics. There has been research on the application of image analysis for phenotyping of an individual plant. Yang et. al adopted x-ray computed tomography to measure rice tiller number. Utilization of digital imaging and image analysis, Granier et al. estimate growth rates from projected leaf area for Arabidopsis plants [START_REF] Christine | an automated platform for reproducible phenotyping of plant responses to soil water deficit in Arabidopsis thaliana permitted the identification of an accession with low sensitivity to soil water deficit [J][END_REF]. Golzarian et. al [START_REF] Mahmood | Accurate inference of shoot biomass from high-throughput images of cereal plants[END_REF]used the projected shoot area of the cereal plants on images at three orthogonal views and plant age in days after planting as predictors to infer the cereal biomass. Jones et al. quantified plant responses to water stress from thermal infrared images of crop canopies [START_REF] Jones Hamlyn | Thermal infrared imaging of crop canopies for the remote diagnosis and quantification of plant responses to water stress in the field[END_REF].

Final yield is a significant agricultural trait in rice breeding. In the screening and evaluation of the rice varieties, measuring and evaluating rice yield is essential. Conventional means of measuring rice yield mainly depend on manual determination. In this method, spikelets were threshed from the panicles and filled spikelets were weighed. Generally, manual measurement is tedious, labor-intensive, subjective and error-prone, especially when large-scale plants were to be investigated.

There has been work on yield-related trait measurement [START_REF] Lingfeng | Fast discrimination and counting of filled/unfilled rice spikelets based on bi-modal imaging[END_REF][START_REF] Lingfeng | A novel machine-vision-based facility for the automatic evaluation of yield-related traits in rice[END_REF]. In these works, rice panicles were harvested from the plant and spikelets were threshed manually or by the threshing machine. The major limit is that, the threshing machine had a threshing error of 3.33% for filled spikelets and 2.27% for total spikelets, and was prone to break the spikelets. Precision agriculture using remote sensing technology, in which the yield of the above-ground canopy for a large area is estimated from satellite and airborne images, is a wide-spread method for yield prediction [START_REF] Yi-Ping | Large-area rice yield forecasting using satellite imageries[END_REF][START_REF] Peijuan | Yield estimation of winter wheat in the North China Plain using the remote-sensing-photosynthesis-yield estimation for crops (RS-P-YEC) model [J][END_REF][START_REF] Jé Go Guillaume | Using Leaf Area Index, retrieved from optical imagery, in the STICS crop model for predicting yield and biomass of field crops[END_REF][START_REF] Anatoly | Use of a green channel in remote sensing of global vegetation from EOS-MODIS[END_REF]. However, this technology is not capable of estimating yield of an individual plant. For yield estimation of individual plant in vivo, academic publication is unavailable.

The objective of the present study is to develop an in vivo, automatic and high-throughput method to estimate the yield of individual pot-grown rice plant using multi-angle RGB imaging and image analysis. In this work, we demonstrated a model that used mixed variables of projected panicle area, projected area of leaf and stem, and fractal dimension achieved good performance in predicting rice yield. In order to improve the generalization ability of our method, we tested our method on rice plants grown both in-door and out-door and achieved promising results.

Experiments and Methods

Rice sample preparation

390 pot-grown rice plants, including 302 rice plants growing in the greenhouse (greenhouse rice) and 88 rice plants growing outdoor (outdoor rice), were used in this study. The samples were harvested at maturing stage and rice yield were measured manually. Yield varied from 7.04 g-43.10 g for the greenhouse rice and 23.68 g-58.43 g for the outdoor rice. The average yield for the greenhouse rice and outdoor rice were 27.55 g and 41.20 g, respectively. The samples were randomly split into two datasets at ratio of 2:1: a training dataset for model construction (260 plants) and a validation dataset for model validation (130 plants).

Hardware setup and image acquisition

Previously, our group developed a high-throughput facility dubbed H-SMART for measuring rice tiller number [START_REF] Wanneng | High-throughput measurement of rice tillers using a conveyor equipped with x-ray computed tomography[END_REF]. The H-SMART facility used an industrial conveyor to transfer pot-grown rice plant to the imaging area for image acquisition. A rotation platform enabled the rice plants to be lifted and turned.

A barcode scanner read barcode of each pot for indexing purpose. To allow measurement of rice yield, we incorporated visible light imaging into the H-SMART facility. As occlusion becomes problematic in rice plants at maturing stage, multi-angle imaging was adopted in the research. Plants were illustrated by fluorescent light tubes both from side and top. Images were taken at every 30°intervals by a Charge Coupled Device (CCD) camera (Stingray F-504C, Applied Vision Technologies, Germany) as the plant was rotated. For each rice plant, 12 images (2452×2056) at different angles were taken.

Lighting conditions were constant throughout the research. Image acquisition was performed by NI-IMAQ Virtual Instruments (VI) Library for LabVIEW (National Instruments Corporation, USA).

More details concerning the H-SMART system can be found in Yang et al. 2011 [START_REF] Wanneng | High-throughput measurement of rice tillers using a conveyor equipped with x-ray computed tomography[END_REF].

Image analysis and feature extraction

Glozarian et al. proved projected shoot area on two dimensional images was a good predictor of shoot biomass of the cereal plants [START_REF] Mahmood | Accurate inference of shoot biomass from high-throughput images of cereal plants[END_REF]. And rice yield was the product of the harvest index and the biomass (more precisely, the total dry weight) [START_REF] Shouichi | Fundamentals of rice crop science[END_REF]. It can be inferred that projected shoot area should have a good correlation with yield. Panicle was the organ where grains grow on. So projected panicle area may correlate well with yield. For this reason, we raised a hypothesis that separating the projected shoot area into two parts: projected panicle area and projected area of leaves and stems would reduce the bias in yield estimation. Gong demonstrated the feasibility of predicting yield per m2 for field rice from fractal dimension and texture features of the rice image [START_REF] Hongju | Estimating paddy yield based on fractal and image texture analysis[END_REF]. In this study, we developed a model using projected panicle area, projected area of leaves and stems, fractal dimension and texture features.

The image analysis and feature extraction was performed in LabVIEW(National Instruments Corporation, USA). Fig. 1 shows the image analysis pipeline for this study. Firstly, the background was removed from the image to get the binary image of rice plant, from which projected area of rice plant (A) and information dimension (IFD) were calculated. Intensity image of the RGB image was extracted according to Eq. 1. Using the binary image as the mask, intensity component of rice plant, from which texture features including differential box counting dimension (DBC) [START_REF] Nirupam | An efficient differential box-counting approach to compute fractal dimension of image[END_REF] and 5 histogram features were calculated, can be obtained. 5 histogram features include the mean value (M), the standard deviation (S), the third moment (mu3), the uniformity (U) and the entropy (E). After that, pixel discrimination and region recognition were performed to obtain the binary image of panicle, from which projected panicle area (PA) was calculated.

I=(R+G+B)/3 (1)

where R, G, B was the R, G, B component of the RGB image. In the end, regions whose area was less than the predefined threshold were removed.

Projected area of rice plant (A) was computed by counting the number of foreground pixels in the binary image of rice plant.

ExG = 2Ng-Nr-Nb (2) i2= (R-B)/2 (3)
where Ng, Nr, Nb was the normalized r,g,b component, defined by Eqns. 4-6.

Nr=R/(R+G+B)

Ng=G/(R+G+B)

Nb=B/(R+G+B) (6)

Fig.2. Image processing diagram for removing background

Fractal dimension can be used to describe shape complexity [START_REF] Bruno Odemir Martinez | Fractal dimension applied to plant identification[END_REF]. The most commonly used method for calculating fractal dimension is box-counting method as it is simple to compute. In this research, the authors adopted information dimension, as it provides a more precise estimate of the fractal dimension than the box-counting method and meanwhile is easy to compute. IFD is calculated by plotting information I(ε) (defined as Eq. 7) against the natural logarithm of box size ε and IFD is computed as the slope of the regression line.

1 1 ln N i i i Ip p      (7)
Where N  is the number of boxes, ε is the box size, i p is the probability of foreground pixels falling into the ith box. Fig. 3 indicates the image processing diagram for generating binary image of panicle. Pixel discrimination between panicle and stem/leaf (non-panicle) based on discriminant analysis was performed to segment the panicle region [START_REF] Hiroshi | Green citrus detection using hyperspectral imaging[END_REF]. There may be non-panicle region in the extracted panicle regions after pixel discrimination. So, region recognition was used to remove non-panicle regions. Rice yield is mainly contributed by effective panicles, which generally locates in the uppermost one third region of the rice plant. Therefore, in this study, a simple location discrimination method was used to recognize the panicle region. More specifically, the image was equally divided into three sub-regions according to the plant height. Regions that were located in the uppermost region were regarded as panicle region and other regions were treated as non-panicle regions. Projected panicle area (PA) was computed by counting the number of foreground pixels in the binary image of panicle. Projected area of leaf and stem (LSA) was computed by A minus PA. Linear regression analysis is one of the most widely used methods in agriculture [START_REF] Mahmood | Accurate inference of shoot biomass from high-throughput images of cereal plants[END_REF][START_REF] Johannes | Evaluation of regression models for above-ground biomass determination in Amazon rainforest [J][END_REF][START_REF] Mf | Allometric models for non-destructive leaf area estimation of Jatropha curcas[END_REF]. In this research, linear regression was used to predict rice yield. Before modeling, all the data were log-transformed in order to meet normality assumptions and increase model performance [START_REF] Cristian | Temperature-controlled organic carbon mineralization in lake sediments[END_REF]. For variable selection and optimal model identification, all-possible regression method was performed, using criteria of Akaike's information (AIC) [START_REF] Hirotugu | A new look at the statistical model identification[END_REF], adjusted coefficient of determination (adjusted R2), and prediction error sum of squares (PRESS statistic) [START_REF] Douglas | Introduction to linear regression analysis[END_REF]. Adjusted R2 indicates the proportion of variation explained by the predictors and increases on adding a variable only when the variable decreases the residual mean square. Taking the complexity of the model into account, AIC is generally regarded as a measure of the goodness of fit of the model. PRESS statistic indicates the prediction ability of the model. A model with a highest adjusted R2 and lowest AIC and PRESS statistic was selected as the final model. Other assumptions requisite to meet the specification of the multiple regression analysis, such as significance of the regression, significance of individual regression coefficients and multi-collinearity, were also checked. Model construction was performed in SAS (version 9.2, SAS Institute Inc., USA).

For model validation, the final model was tested on the validation dataset. In addition, cross validation technique was used to evaluate the prediction error of the model [START_REF] Mahmood | Accurate inference of shoot biomass from high-throughput images of cereal plants[END_REF]. In the cross validation method, the samples are randomly partitioned into K approximately equal-sized K groups. In each iteration, the i-th (i=1,2…K) group is treated as validation set, and the remaining K-1 groups are used to fit the model. This process are repeated for K iterations, and the prediction errors, such as root mean square error (RMSE) and mean absolute percentage error (MAPE), are averaged over the iterations. Model validation was performed in Matlab ((Mathworks Inc., USA).

Results and Discussion

Image analysis

Rice yield is mainly contributed by effective panicles, which generally locates in the uppermost one third region of the rice plant. In this study, the candidate panicle regions were firstly segmented using discriminant analysis. Then, the panicle regions were discriminated from non-panicle regions according to their location. After that, projected panicle area (PA) was computed. This method would be fit for most rice varieties. However, some rice varieties may have effective panicle that locates in the middle of the plant. In this occasion, some panicle regions may be mistaken treated as non-panicle regions and consequently the computed projected panicle area may be biased. This paper mainly focuses on demonstrating the feasibility of estimating rice yield from projected panicle area, projected area of leaf and stem and fractal dimension. Improvement of projected panicle area extraction would be our future orientation.

Model development and validation

In variable selection and optimal model identification, we noticed that when the number of variables in the model was fixed, the model with the highest R2 generally had the lowest AIC and PRESS statistic.

Therefore, among models with the same number of variables, the model with the highest R2 was selected for further consideration to accelerate the optimal model identification. This process reduced the number of candidate models to nine. Table 1 shows the performance of the nine candidate models.

When using alone, IFD was the most powerful variable for predicting yield, explaining 62.8% of the variation in the 260 training samples. Combination of PA and LSA outperformed other models with two variables, increasing adjusted R2 to 0.87, and dramatically reducing AIC and PRESS. When the number of variables in the model was more than 4 (excluding constant), adding more variables to the model did not produce less AIC and PRESS. The model with four variables IFD, E, PA and LSA had the highest Adjusted R2 and lowest AIC and PRESS. However, multi-collinearity was noticed in this model (VIF>5).Therefore, the model with three variables (Eq.8 and Table 2) was selected as the final model. Model diagnostics suggested that underlying assumptions of the regression model were met. 2 shows the summary of the final model based on the training samples. F-test value and its significance (P>F) showed that the overall fit of the model was significant. The t-test values and its significance (P>|t|) indicated that the individual independent variables were statistically significant.

Model 1: ln(Y)=a 0 + a 1 • ln(PA)+a 2 • ln(LSA)+a 3 • ln(IFD) (8)
Variance inflation factor (VIF) for each individual independent variable showed that there was no multi-collinearity between the variables. The final model with three variables was then applied to the validation dataset and performed well (Fig. 5). The RMSE and MAPE for the validation dataset were 3.04 g and 7.74%.

5-fold cross validation showed that the RMSE and MAPE of the final model were 2.83 g and 7.45%.

Performance comparison of the different models

For comparison, we also tested the model using projected shoot area (Eq.9) and the model using projected panicle area and projected area of leaf and stem (Eq.10). And compared the two models with the model used in this study (Table 3). Separating A into PA and LSA significantly reduces the estimation error in yield prediction, which proved our hypothesis. And Adding fractal dimension into the model had a slight promotion on model performance. This was because that fractal dimension describes shape complexity of rice and can compensate the influence of organ occlusion problem at a certain degree.

Performance comparison of the models using different numbers of side-view images

By increasing number of side-view images, more hidden organs were available to be seen and more complete information on plant architecture would be provided. To find out the appropriate number of side-view images, the models that used average prediction variables from 1, 2, 4 and 6 side-view images were studied. The four models were listed as follows. Table 4 shows the performance comparison of the models using different numbers of side-view images.

5-fold cross validation indicated that increasing number of side-view images would decrease prediction error of the model. However, the prediction error reduction was not significant when the number of side-view images were larger than 4. In this study, the images were captured by the visible light imaging unit that incorporated into our previous H-SMART system. The processing time of the H-SMART system was 20 seconds. This time was sufficient for capturing and processing 12 side-view images. For this reason, 12 side-view images were used in this study. To promote the efficiency of the yield estimation, two perpendicular side-view images would generate acceptable results.

Conclusions

This paper presented a non-destructive and high-throughput method to estimate the yield of individual pot-grown rice plant using photonics-based technology. And we demonstrated a new idea of estimating rice yield from projected panicle area, projected area of leaf and stem and fractal dimension. 5-fold cross validation showed that the predictive error was 7.45% with testing rice plants grown both in-door and out-door. The presented work has the potential of accelerating yield estimation and would be a promising impetus for plant phenomics and high-throughput phenotyping. (5) Where G i was the i-th graylevel, and p(G i ) was the probability of G i .

Acknowledgment

Fig. 1 .

 1 Fig.1. Image analysis pipeline for this study Fig.2 illustrates the image processing diagram for removing background. First, the ExG component and i2 component were extracted according to Eqns. 2-3. The ExG component was used to segment the green parts of the rice plant, and the i2 component was used to segment the yellow parts. The images inthis study were captured in a constant illustration environment, so both the ExG component and i2 component were binarized using fixed thresholds. Then the union of the two binary images was calculated. In the end, regions whose area was less than the predefined threshold were removed.

Fig. 3 .

 3 Fig.3. Image processing diagram for removing background

Fig. 4

 4 Fig.4 shows the plot of residuals versus fitted values ln( ) Y of the final model. It indicated that the residuals were evenly distributed at the two sides of y=0 and can be contained in a horizontal band,suggesting no obvious model defects were existed. In addition, no significant differences were presented among outdoor rice and greenhouse rice. This meant that the model performed well on both outdoor rice and greenhouse rice.

Fig. 4 .

 4 Fig.4. Plot of residuals against the fitted values ln( ) Y of the final model

Fig. 5 .

 5 Fig.5. Prediction for testing set using the final model

Model 2 :

 2 ln(Y)=a 0 +a 1 • ln(A)[START_REF] Lingfeng | Fast discrimination and counting of filled/unfilled rice spikelets based on bi-modal imaging[END_REF] Model 3: ln(Y)=a 0 +a 1 • ln(PA)+a 2 • ln(LSA)[START_REF] Lingfeng | A novel machine-vision-based facility for the automatic evaluation of yield-related traits in rice[END_REF]

Model 4 :

 4 a single image at angle of 0 Model 5: two images at the angle of 0 and 90 Model 6: four images at the angle of 0, 90, 180 and 270 Model 7: six images at the angle of 0, 60, 120, 180, 240 and 300

  This work was supported by grants from the Program for the National Program on High Technology Development (2013AA102403), the National Natural Science Foundation of China (30921091, 31200274), and the Fundamental Research Funds for the Central Universities (2013PY034), New Century Excellent Talents in University (No. NCET-10-0386). Additional equations: Definition of the five histogram features.

Table 1 .

 1 Performance of the nine candidate model based on training samples

	Number of variables Adjusted R 2 AIC	PRESS Variables in the model
	1	0.6284	-940.80	7.01	IFD
	2	0.8725	-1218.04 2.40	PA, LSA
	3	0.8904	-1256.28 2.08	IFD, PA, LSA
	4	0.8917	-1258.37 2.06	IFD, E, PA, LSA
	5	0.8916	-1257.08 2.07	DBC, IFD, E, PA, LSA
	6	0.8917	-1256.44 2.08	IFD, M, S, mu3, PA, LSA
	7	0.8915	-1255.06 2.10	IFD, E, M, S, mu3, PA, LSA
	8	0.8915	-1253.94 2.11	IFD, E, M, S, U, mu3, PA, LSA
	9	0.8910	-1251.98 2.13	DBC, IFD, E, M, S, U, mu3, PA, LSA
	Table

Table 2 .

 2 Summary of the final model based on training samples

	Coefficients Value Standard error	t	P>|t|	VIF	F	P>F	R 2	RMSE
	a 0 -10.613	0.39223	-27.06 <.0001	0
	a 1 0.4181 a 2 0.64369	0.01892 0.0439	22.1 14.66 <.0001 2.09423 <.0001 1.36293	702.32 <.0001 0.89 2.71 g
	a 3 3.67453	0.56133	6.55	<.0001 2.59348

Table 3 .

 3 Performance comparison of the three models

	Training set		Validation set	Cross validation
	RMSE (g) MAPE RMSE (g) MAPE RMSE (g) MAPE
	Model 1 2.71	7.16%	3.04	7.74%	2.83	7.45%
	Model 2 4.79	12.03% 5.48	13.82% 5.04	12.74%
	Model 3 2.94	7.76%	3.35	8.48%	3.09	8.08%

Table 4 .

 4 Performance comparison of the models using different numbers of side-view images

	Model	Number of side-view images	5-fold Cross validation RMSE (g) MAPE
	Model 4 1	3.80	10.08%
	Model 5 2	3.13	8.15%
	Model 6 4	3.01	7.87 %
	Model 7 6	2.92	7.77%
	Model 1 12	2.83	7.45%


Jing Wang

Kai Xia

Xuegong Chen

Research and Implementation of modeling Grid DEM based on discrete data

Keywords: DEM interpolation, digital elevation model, distance-weighted average method, least squares fitting method

Interpolation method is used to construct regular grid DEM (Digital Elevation Model) based on discrete data. The key to DEM interpolation algorithm lies in the efficiency of data search around the points to be operated and the precision of DEM interpolation by selecting an appropriate interpolation function model. At present, the point by point interpolation is very common. Based on the comparison and analysis of characteristics and application scope of various DEM interpolation algorithms, this paper proposed a new interpolation model which combines distance-weighted average method and least squares fitting method. Results show that proposed method can improve search speed and interpolation accuracy remarkably and has good application value for large area.

Introduction

With the rapid development of GIS (Geographic Information System), DEM (Digital Elevation Model), as a digital description and simulation of the surface topography, has become an important component of spatial data infrastructure and "Digital Earth". The accuracy of DEM has a significant impact on terrain visualization, GIS analysis and credibility of the decision-making [START_REF] Zhilin | Digital Elevation Model (Second Edition)[END_REF][START_REF] Roger H Twito | The MAP Program:Building the Digital Terrain Model[END_REF]. It is necessary to improve the accuracy of modeling DEM. Most of original data are discrete data. The key issue of DEM is how to construct regular grid DEM by interpolating with discrete data. Therefore, research on improving the DEM accuracy is of great significance.

To construct Grid DEM, first we divide the research region into two-dimensional grid to form spatial structure of grid that covers the entire area. Then discrete points surrounding the interpolated grid points are used to calculate its elevation values, we can obtain the grid data of this region with a certain format. Interpolation is the core issue of modeling DEM, which involved in DEM in all aspects of production, including quality control, accuracy assessment and analysis applications. According to the distributed range of interpolation points, DEM data interpolation algorithm generally can be divided into overall interpolation, block interpolation and point by point interpolation. The point by point interpolation is widely used because its calculation method has many advantages, such as simple, flexible block, high precision and so on. The point by point interpolation can be divided into the distance-weighted average method, moving surface fitting method, least squares fitting method, etc. Every interpolation method has some drawbacks. For example, the distanceweighted average method is less precise and poor in smoothness, and only suitable for points which are well-distributed and sampled from even terrain; Moving surface fitting method is of high precision and flexible in calculation, but imposes high requirements on the sampling point; Least squares fitting method is suitable for most of terrain but time consuming [START_REF] Tang Guoan | Digital Elevation Model (Second Edition)[END_REF][START_REF] Charleux Demargne | Quality assessment for drainage networks and watershed boundaries extraction from a digital elevation model (DEM)[END_REF][START_REF] Aguilar | Effects of terrain morphology, sampling density, and interpolation methods on grid DEM accuracy[END_REF][START_REF] Gl Heritage | IC. Fuller influence of survey strategy and interpolation model on DEM quality[END_REF][START_REF] Lazzaro | Radial basis functions for the multivariate interpolation of large scattered data sets [J][END_REF][START_REF] Li | A simple efficient algorithm for interpolation between different grids in both 2D and 3D[END_REF]. Different interpolation algorithms have their own advantages, disadvantages and application scope, so we should select the appropriate interpolation algorithm according to a certain problem.

Simply using one method of interpolation to construct DEM in entire region have limitations, such as high requirements on the sampling point, time consuming and so on, which will affect the applicability and accuracy of model. People put forward some combination method, for example, integration of moving surface fitting method with distance-weighted average method, distance-weighted average method with quadratic surface method. Both of these methods only suitable for points which are well-distributed and sampled from even terrain. The method that combines moving surface fitting method with least squares fitting method does not consider moving surface method has a high requirement on accuracy of sampling data [START_REF] Pan | Study and practice of DEM generation [J][END_REF][START_REF] Yin | Improvement of DEM interpolation algorithms based on moving surface fitting and distance-weighted[END_REF][START_REF] Jianming | Integration of the two methods to interpolate grid DEM [J][END_REF][START_REF] Hu | The study of DEM modeling method based on regular grid[END_REF][START_REF] Jinming | Effects of interpolation parameters in inverse distance weighted method on DEM accuracy[END_REF].

In response to above problems, taking all conditions into account and after analysis of the application of various interpolation methods, this paper presents a method that integrates the distance-weighted average method with the least-squares fitting to solve problems arising in DEM modeling and improve accuracy and efficiency of the interpolation method.

The following of the paper is organized as follows. Section 2 presents the proposed approach. Section 3 describes an experiment and shows the analysis. Section 4 concludes the paper.

2

Data and Methods

Data

These elevation data used in this paper are sampled from three areas, such as alpine areas, glacier areas and even areas. Each area contains more than 50000 sampling points. Using the distance-weighted average method, least squares method, moving surface fitting method, moving surface fitting method and the proposed combination method to construct DEM based on discrete data, we analyze and evaluate the accuracy of the model. Set aside 100 sampling points as the checkpoints that do not participate in DEM modeling. Using different methods to model DEM, comparing the checkpoints original elevation values and calculated elevation values through four interpolation algorithm to calculate max error and RMSE (root mean square error), then we analyze and evaluate the results.

Distance-weighted average method

To use the distance-weighted average method, we need to select these discrete points nearest to the interpolation point. The easiest way is to directly calculate the distance between interpolation points and the discrete points, then choose discrete points according to distance and interpolate, but this interpolation method need large time and space to compute, which has a bad impact on interpolation speed [START_REF] Yin | Improvement of DEM interpolation algorithms based on moving surface fitting and distance-weighted[END_REF]. So we choose certain search areas according to requirement of mathematical model for the interpolation number and adjust the search scope.

The distance-weighted average method is as follows:

     m k k m k k k w h w 1 1 h (1)
h k is the point to be interpolated, w k is the weight of discrete point k, and h k is elevation values of sampling point k.

Currently weight calculated in the following ways:

-u ki i d P  (2) 2 i P           ki ki d d R (3) 2 i e P        m d ki (4)
R is a search radius; m is a constant, P i is the weight of discrete points. In this paper, we choose the first method, and set u = 2.

A problem arises when there are few reference points. When the reference points are very far from each other or distributed unevenly, if we still use these reference points that has little effect on the interpolated points to interpolate, it will cause disturbance of interpolation points instead of improve the accuracy of the interpolation points, which make the interpolation results not ideal. To improve the efficiency of distance-weighted average algorithm, this paper adopts the idea of gradually expanding the search area and restricted orientation search. The idea is as follows:

First, classify grid in accordance with the actual needs. Second, determine the number of discrete points needed to be calculated.

Third, search discrete points around the interpolated grid points with the restricted orientation. If there are enough discrete points in each quadrant, use the distanceweighted average method for interpolation. If not, expand search area gradually. When search area reaches threshold, then stop.

Least squares fitting method

When the reference points are very far from the interpolation points or distributed unevenly, distance-weighted average method has great error. Least squares fitting method, which combines the advantage of distance-weighted average method and trend surface method, is suitable for most of terrain, and can effectively eliminate noise data [START_REF] Hu | The study of DEM modeling method based on regular grid[END_REF]. Weighted least squares fitting method assigns all discrete points that have influence on interpolation points a certain influence coefficients, and assumes the coefficients is consistent with a quadratic polynomial as follows:

2 6 2 5 4 3 2 1 c c y) p(x, y c x c xy c y c x       (5)         2 2 n 2 i i ] , x [p Q b y a x w z y i i i i        (6)
When we use least squares fitting method for interpolation, these points that are nearer to the interpolation point are allocated a higher weight than these points which are far away. As the number of discrete points the least squares fitting method requires is more than 6,we need to find coefficients ci ( i= 1, 2, … , 6 ) to make (6)

minimum.       2 2 b y a x w i i    is weight and is a function of the distance. Take             2 2 2 2 2 1 b y y x b y a x w i i i i i       
, when (a, b)close to (x i , y i ), it has a large value and a small value while (a, b)is far from (x i , y i ). To get minimal Q, according to the principle of least squares

0 i c Q   
, so we can get the following equation:

    0 , 2 c Q 1 i 1        w z y x P N i i i (7)     0 , 2 c Q 1 i 2        i N i i i wx z y x P (8)     0 , 2 c Q 1 i 3        i N i i i wy z y x P (9)     0 , 2 c Q 1 i 4        i i N i i i y wx z y x P (10)     0 , 2 c Q 2 1 i 5        i N i i i wx z y x P (11)     0 , 2 c Q 2 1 i 6        i N i i i wy z y x P (12)
The equations from ( 7) to ( 12) can be solved to obtain coefficients of p (x i , y i ). Use these coefficients in equation ( 5); we can obtain value of grid points. Solving these equations needs large time. When there are large discrete points, there is a sharp decline in speed. When the reference points are very far from the interpolation points or distributed unevenly, least squares fitting method has better accuracy than distanceweighted average method. So we combine least squares fitting method and distanceweighted average method to improve the accuracy and efficiency of interpolation.

Proposed method

When the reference points are very far from each other or distributed unevenly, the result of using weighted average algorithm would not be ideal. This paper brings up the idea of gradually expanding the search area and restricted orientation search. As the number of discrete data points the distance-weighted average method requires is from 4 to 8, we adopt four quadrants search method, so there should be at least one point in each quadrant. The least squares fitting method has no special requirements on the distribution of discrete points, but it needs at least six discrete points. When the number of discrete points is over six and the distribution does not satisfy the requirements of distance-weighted average method, we adopt least squares fitting method.

In addition, in order to accelerate the computing speed and improve the efficiency of modeling, we judge distance between point to be interpolated and discrete point firstly. For example, if the nearest distance between point to be interpolated and discrete point is less than 0.5 meter, then we use value of the elevation of discrete point directly and don't do next step, if not, use the proposed algorithm for interpolation. Algorithm is as follows:

First, for point P to be interpolated, search discrete points in grid with number 1 around the point to be operated with the restricted orientation searching. If each of four quadrants has more than one discrete point, use distance-weighted average method for interpolation; Second, if not, expand search area with number 2. Check if there is more than one discrete point in each quadrant, if yes; use distance-weighted average method for interpolation; Third, if discrete points in one or more quadrants are less than one, but a total of discrete points are more than six, use least squares fitting method to interpolate.

Fourth, if discrete points are less than 6, expand the search area with number 3 and use least squares fitting method to interpolate if the requirement is satisfied ; Sixth, if not, don't do interpolation and process the next grid point.

Experiments

Using distance-weighted average method, least squares fitting method, moving surface fitting method and proposed method to model DEM. Set aside 100 sampling points as the checkpoints that do not participate in DEM modeling, if the check point distribute unevenly, we may get an inaccurate evaluation of interpolations, in this paper we make the check point distribute evenly in order to obtain an accurate evaluation. We compare the original elevation values and calculated elevation values through four interpolation algorithms to calculate RMSE and max error. In order to directly display the results, we randomly selected 10 sampling points from 100 sampling points to compare. Fig. 2, Fig. 3 and Fig. 4 show absolute value of error in alpine areas, glacier areas and even areas respectively. Table 1, table 2 and table 3 show comparison of various interpolation precision in alpine areas glacier areas and even areas respectively. The results are as following:

Discussion

From these tables, we can sum up:

The proposed method has higher accuracy than distance-weighted average method and least squares fitting method:

In alpine areas, RMSE of using proposed method is 8.48, RMSE of using least squares fitting method is 12.06, RMSE of using distance-weighted average method is 15.27, and RMSE of using moving surface fitting method is 14.11, an decrease of 3.58, 6.79, 5.63 respectively.

In glacier area, RMSE of using proposed method is 6.54, RMSE of using least squares fitting method is 9.91, RMSE of using distance-weighted average method is 11.36, and RMSE of using moving surface fitting method is 12.24, an decrease of 3.37, 4.82, 5.7 respectively.

In even areas, RMSE of using proposed method is 3.37, RMSE of using least squares fitting method is 4.46, RMSE of using distance-weighted average method is 5.78, and RMSE of using moving surface fitting method is 4.51, an decrease of 1.09, 2.41, 1.14 respectively.

This result amply demonstrated that the proposed method has significantly improved accuracy in various areas, especially in alpine areas. This proved that the proposed method has a great of advantage in meeting actual demand and particularly suitable for the modeling of uneven terrain of the area.

The distance-weighted average method and moving surface fitting method has weak accuracy in alpine areas and glacier areas, better accuracy in even areas. The least squares fitting method has higher accuracy than distance-weighted average method in each kind of areas, especially in undulating areas. But in even area, accuracy is only a few higher. Taking the computing time into account, the least squares method is suitable for terrain that changes significantly.

Conclusions

Using distance-weighted average method, least squares fitting method, moving surface fitting method and proposed method to model DEM, according to the experimental results, after analysis, we can get the following conclusions.

The proposed method which combines distance-weighted average method and least squares fitting method has adopt both advantages and improved accuracy in each kind of terrain; it is suitable for the modeling of large areas. The experiment has demonstrated that the proposed method has improved efficiency of interpolation and prevented generating unwanted interpolation points.

The thinking of distance-weighted average method takes large ground gradient as a continuous process. But when the terrain changes significantly, such as ridges or rifts, the accuracy is not particularly ideal. So this method is suitable for even areas.

The least squares fitting method has high accuracy and time consuming, so it is suitable for terrain changes significantly.

Moving surface fitting method is suitable for points which are well-distributed and sampled from even terrain. When the reference points are very far from each other or distributed unevenly, this method would not be ideal, because it has a high requirement on accuracy of sampling data.

It is difficult to find a model that fully fit the complex topography of the surface. Therefore, in practical applications, how to make better use of space to fit the data and actual ground is an important direction for future research.

Fig. 1 .

 1 Fig. 1. Grid model searching

Fig. 2 .Fig. 3 .Fig. 4 .

 234 Fig. 2. Absolute value of error in alpine areas

Table 1 .

 1 Comparison of various interpolation precision in alpine areas (unit: meter)

	Interpolation method	Distance-weighted average method	Least fitting method squares	Moving fitting method	surface	Proposed method
	RMSE	15.27	12.06	14.11		8.48
	Max error	28.67	20.64	22.53		13.77

Table 2 .

 2 Comparison of various interpolation precision in glacier areas (unit: meter)

	Interpolation method	Distance-weighted average method	Least fitting method squares	Moving fitting method	surface	Proposed method
	RMSE	11.36	9.91	12.24		6.54
	Max error	23.13	18.71	27.15		12.92

Table 3 .

 3 Comparison of various interpolation precision in even areas (unit: meter)

	Interpolation method	Distance-weighted average method	Least fitting method squares	Moving fitting method	surface	Proposed method
	RMSE	5.78	4.46	4.51		3.37
	Max error	10.81	9.25	9.33		8.16


Zilong Wang

Qiang Fu

Jun Meng

Qiuxiang Jiang

Xianghao Wang

Research on Spatial Variability Characteristics of Black Soil Unfrozen Water in Songnen Plain during Freezing-Thawing Period

Keywords: Songnen Plain, black soil, soil unfrozen water, spatial variability

To address the difficulties of the complexity and the quantitative description of spatial non-homogeneity in the process of freezing and thawing of soil moisture migration, taking the unfrozen water as the research object, and using geostatistics to study the spatial variability characteristics of unfrozen soil water in different freezing and thawing stages. The results show that different freezing and thawing stages of unfrozen soil water has good spatial structure and strong spatial correlation, and it has a strong spatial redistribution effect on soil moisture. The application of geostatistics provides a new train of thought for in-depth study on freezing-thawing process of soil moisture transport mechanism.

Introduction

Freezing and thawing of soil water movement, as an important part of the water cycle in nature, agriculture, water resources, environmental engineering, which occupies an extremely important position. To carry on an in-depth investigation into the variations of the freezing and thawing period soil moisture, master the laws of its motion, will not only help to promote the development of the theory of soil water dynamics, but also provide the theoretical basis for water resources accurately evaluation, efficient utilization of soil water, and a reasonable determination of farmland irrigation technical parameters, etc [START_REF] Guodong | Recent Development of Geocryological Study in China [J][END_REF].

Study on the theory of soil water dynamics have made considerable progress, along with the changes of soil moisture from the morphology point of view to the energy state theory [START_REF] Zhidong | Soil water dynamics [M[END_REF],the research methods also vary from the qualitative description [START_REF] Bloomsberg | Effect of Moisture Content on Permeability of Frozen Soils [C[END_REF][START_REF] Guoqing | Freezing Point of Loess in Lanzhou, China[J][END_REF] based on soil physics test to the quantitative research [START_REF] Shoop | Moisture Migration during Freeze and Thaw of Unsaturated Soils: Modeling and Large Scale Experiments[END_REF][START_REF] Xiuqing | Numerical Emulation on Simultaneous Soil Moisture and Heat Transfer under Freezing and Thawing Conditions[END_REF] based on the numerical simulation. However, dramatic phase transitions and spatial variability of soil moisture in the freezing thawing process, that lead to the existing methods in the analysis of soil moisture variation during freezing and thawing period insufficient.

Geostatistics is a theory, which study on the spatial variability of the properties of a system. The theory has been applied in the study on spatial variability of in unfreezing soil water [START_REF] Qiuxiang | Research on Spatial Variability of Soil Water Characteristics In Western Semi-arid Area of Heilongjiang Province [J][END_REF].But it has fewer applications on the freezing and thawing soil moisture variation.

Therefore, in this paper, the methods of geostatistics were used to analysis the spatial variability of black unfrozen soil water in Songnen Plain in the freezingthawing period, expecting to provide scientific basis for reasonable regulation and utilization of black soil water in Songnen Plain.

2

Materials and Methods

Study Area and Test Scheme

Field experiments were measured in a place of dry land in the city of Harbin, soil type is black soil. Harbin city is located in the southeast of Songnen Plain, has a temperate continental monsoon climate with its severe cold winter, dry and windy spring, has no perennial frozen soil layer, the annual average temperature is 3.6℃, the average January temperature is -20.3℃, the average frost free period is 141 days, the average annual precipitation is 523.3mm, the average annual snowfall is 63.1mm, the annual average frost depth is about 1.75m. Field experiment treats the soil unfrozen water in freezing-thawing stages as the observation object, and uses TRIME-T3 time domain reflectometry to monitor soil unfrozen water during the soil thawing period of 2008 to 2009. In the same vertical profile every 2m buried 1 TDR probe, altogether 10 probes was buried, each TDR probe sets 8 measuring points, respectively 20cm, 40cm, 60cm, 80cm, 100cm, 140cm, 180cm, 220cm and 260cm.

Research Methods

Geostatistics is the most important method to study the spatial variability characteristics of soil. Since 1963, the famous French statistician Matheron G has put forward the concept of geostatistics [START_REF] Matheron | Principles of Geostatistics[END_REF], through decades of development, geostatistics has been widely applied in soil science, ecology, geology, hydrology, meteorology, environment and many other fields.

Since geostatistics methods and theory can not only keep the autocorrelation of space variables, but also can make the simulation value of measured point be equal to the measured value, it can use a few sampling point to simulate many unknown information and spatial distribution, which greatly reduced the field sampling and observed workload, enhanced visibility [START_REF] Yaxin | Forecast Theory and Conditional Simulation of Soil Water and Salt Information's Spatial Variability[END_REF] of data analysis.

3

Results and Analysis

Semivariance Analysis

Based on the soil unfrozen water content observational data in different freezingthawing stages, using the software of GS+ 7 for the semi variance analysis, as shown in Table 1. Results of the Semi variance analysis of soil unfrozen water content in different freezing-thawing stages in the experimental areas show that: soil unfrozen water content in different freezing-thawing stages has good spatial structure, their theoretical semivariogram models were all spherical model, as the same as the theoretical semivariogram models of unfrozen soil total water content in different freezing-thawing stages in previous research [START_REF] Zilong | Spatial Variability of Soil Moisture Profile in Seasonal Frozen Soil Region in Different Stages[END_REF].The spatial structure of the soil moisture test area should adopt the spherical model to describe. Model fitting precision is lower than the total soil moisture content's, but still can meet the precision requirement. As is shown from the unfrozen soil water content Nugget/Sill in different freezing-thawing stages, the spatial correlation of most unfrozen soil water content is stronger than the corresponding position of the unfrozen soil water content [START_REF] Cambardella | Field-scale Variability of Soil Properties in Central Iowa Soils[END_REF].

According to the horizontal time span comparison, the spatial correlation of the initial freezing period (December 19, 2008) and the melt season (April 23 and May 19, 2009) is relatively weak. Soil has strong spatial correlation in stable freezing period. Effect of freezing process on soil unfrozen water content spatial correlation was in contrast to the situation in soil water content.

Investigation of the range of soil unfrozen water content: range of soil unfrozen water content in every freezing stages of the freezing period (before April 23, 2009) was about 14m,the range of melting period decreased by 50%.It shows that the freezing-thawing processes also have impact on the spatial correlation distance of unfrozen soil water content, but relatively small.

Analysis of Spatial Variation Characteristics

Using the Geostatistical Analyst module in ArcGIS 9 Software to carry on the spatial local estimation of the unfrozen soil water content in different freezing-thawing stages and to map the spatial distribution, as shown in Fig. 1.The left side data in the spatial distribution of the unfrozen soil water content in different stages in Fig. 1 represents the depth of the measuring points, the unit is cm; the upper marks represents each TDR probe. From Fig. 1 we can see, in the process of freezing-thawing stages, the two dimensional spatial distribution of the unfrozen soil water content in the profile has certain regularity. The unfrozen soil water content and soil the negative temperature has an extremely close relationship [START_REF] Xuezu | Experimental Study on Water Migration in Frozen Soil[END_REF][START_REF] Zhou Youwu | Chinese Frozen Soil[END_REF]. In December 19, 2008, the soil belongs to the initial freezing period, freezing depth is 65cm, the equivalent zone of soil unfrozen water content present banded rules, it shows that the phase change process of water plays an important role in the redistribution of soil moisture. For the part above 60cm measuring point in the freezing layer and the part below 100cm measuring point in the unfrozen layer, the banded equivalent zone of soil unfrozen water content is wide, that is to say soil unfrozen water in these region changes little in the vertical direction. However, the equivalent zone of soil unfrozen water content at the measuring points between 60cm and 100cm is narrow, it changes heavily in the vertical direction. Because of the existence of soil freezing front, the unfrozen soil water below the location of the freezing front migrate to the frozen zone, this phenomenon dovetails the fact of the freezing depth of 65cm.

Based on the same principle, the freezing depth is about 100cm in January 10, 2009, 160cm in March 16, 2009, and in March 27, 2009 the freezing depth reached the maximum 162cm, in the adjacent soil layer below measuring points for 100cm and the soil layer between 140cm to 180cm between the unfrozen water content dramatic changes in the vertical direction, which were associated with the development of freezing depth or the depth of the freezing front, that's to say that the drastic change in the vertical direction of soil unfrozen water content in each frozen stages was caused by soil freezing. Therefore, the location of the freezing soil can be inferred according to the spatial distribution chart of the soil unfrozen water content in the freezing period. The unfrozen soil water content in freezing period showed a trend of increasing gradually from top to bottom. At the same time, notably in the freezing period of soil, for the frozen layer, the spatial distribution of the unfrozen water content presents regularly banded structure as the same as in the initial freezing period. It shows that when the soil is frozen, the unfrozen water content decreases rapidly, and continues to decline with the temperature's further reducing, but the freezing effect on the movement of the unfrozen water in frozen layer is small, the possible reasons to this question mainly are the atmospheric environmental conditions, the boundary of micro topography and the spatial variability of the soil water conductivity characteristics.

For the soil thawing period, the soil began to melt from two directions [14-15], while the soil unfrozen water content also increased from two directions. The spatial distribution situation of the soil unfrozen water content in the whole soil profile in April 23 and May 19, 2009 were shown in Fig. 1. Frozen soil region is the low value area of unfrozen water content in soil, such as the measuring point between 40cm to 140cm in April 23, 2009 and the measuring point between 100cm to 140cm in May 19, 2009. At the same time, according to the spatial distribution of soil unfrozen water content in soil thawing period melting position between the upper and lower soil layer could be roughly inferred. For the soil surface in melting period, the initial melting period was the main driving force for migration of the unfrozen soil water content. Thus, the distribution of the soil surface unfrozen water content was regularly banded, and so was the situation in initial soil melting period. The soil unfrozen water content showed strong spatial variability after the surface soil began to thaw for a period of time. This is because in this period soil unfrozen water content is mainly affected by the different position of soil evaporation and the spatial variability of the water conductivity characteristics [START_REF] Shi Chengxi | Agricultural Hydrology[END_REF]. The lower soil also melted, but was less affected by external factors, so the unfrozen water content showed no spatial variability characteristics.

In order to analyze the severity that the different soil unfrozen water content vary with time in soil freezing-thawing period, spatial distribution maps of soil unfrozen water content coefficient of variation were drawn, as shown in Fig. 2. We can see from Fig. 2, during the whole freezing-thawing period, coefficient of variation of the upper layer of the soil unfrozen water content is large, the lower layer of the soil unfrozen water content on the contrary. This is due to the upper layer of the soil unfrozen water content is not only the influence of the freezing-thawing effect, but also restricted by other factors. The coefficient of variation of bottom soil unfrozen water content is very small, below 5.57%. The phase change of the soil moisture in freezing period did not occur this part of soil layer, so there were many other factors that bring about the slight variation of its soil unfrozen water content with time changes.

Conclusions

Different stages of freezing and thawing soil not frozen water content has good spatial structure and strong spatial correlation, but has relatively small influence on the spatial correlation distance of the soil unfrozen water content. In the whole process of freezing-thawing, the severe variation of soil water content in vertical direction happened in the position of soil freezing front. The spatial distribution of soil unfrozen water content was mainly influenced by the effect of freezing-thawing in the initial freezing-thawing period, and was influenced by the environment and the variability of the soil itself in the late period of freezing-thawing period.

Using geostatistical theory and method is a good way to analyze different stages of freezing-thawing soil spatial variation features of unfrozen water content, also can visually depict the spatial distribution state's variation of soil unfrozen water content in different freezing-thawing period, at the same time, it can infer the freezingthawing position and the variability of the whole soil profile according to the spatial distribution of the unfrozen soil water content in the freezing-thawing period, so as to provides a new thought and method for further research on soil moisture migration mechanism of Songnen plain's black soil region at the freezing-thawing period.

Fig. 1 .

 1 Fig.1. Spatial distribution maps of soil unfrozen water content in different freezing-thawing stages

Fig. 2 .

 2 Fig.2. Spatial distribution map of coefficients of variation of soil unfrozen water content in freezing-thawing period

Table 1 .

 1 Semivariance analysis of soil unfrozen water content in different freezing-thawing stages

	Date	Fitting model	Nugget	Sill	Nugget/Sill (%)	Range (m)	R 2
	08-12-19	Sphere	0.228	1.323	17.2	14.41	0.827
	09-01-10	Sphere	0.001	1.443	0.1	14.8	0.899
	09-03-16	Sphere	0.001	1.41	0.1	14.2	0.876
	09-03-27	Sphere	0.001	1.392	0.1	13.91	0.873
	09-04-23	Sphere	0.216	1.081	20.0	9.82	0.737
	09-05-19	Sphere	0.238	1.078	22.1	6.47	0.912

Acknowledgment

The study was supported by the National Natural Science Foundation of China (51209039, 51279031), Specialized Research Fund for the Doctoral Program of Higher Education (20122325120016), China Postdoctoral Science Foundation (2012T50321, 2011M500632), Program for New Century Excellent Talents in University, Science Foundation for Distinguished Young Scholars of Heilongjiang Province (JC201402), Natural Science Foundation of Heilongjiang Province (E201241), Foundation of Heilongjiang Educational Committee (12521016), Heilongjiang Postdoctoral Fund (LBH-Z11231) and Doctoral Foundation of Northeast Agricultural University (2012RCB13).


Qiuxiang Jiang
email: jiangqiuxiang914@163.com

Qiang Fu

Jun Meng

Zilong Wang

Ke Zhao

Comprehensive Evaluation of Land Resources Carrying Capacity under Different Scales Based on RAGA-PPC

Keywords: land resources, carrying capacity evaluation, projection pursuit evaluation, Sanjiang Plain

Land carrying capacity is an important indicator of the regional population, resources and environment evaluation of sustainable development. In this study, based on fully considering the influencing factors of land resources carrying capacity, land resources carrying capacity evaluation indicator system is constructed, composed of four sub target, i.e. the utilization of land resources, social development, economic and technology and the ecological environment. Real-coded accelerating genetic algorithm is used to optimize projection pursuit model to carry on comprehensive evaluation of land resources carrying capacity under different scales. Results show that national resources carrying capacity level of large-scale land falls under grade III under national scale; Heilongjiang Province and Sanjiang Plain belong to grade II under medium scale; The city of Jixi, Hegang, Shuangyashan and Jiamusi fall on grade I and the city of Qitai River, Muling and Yilan belong to grade II under small scale. The measures of effective utilization of land resources can be further proposed combined with the regional land resource evaluation results.

Introduction

Land is the material basis for human survival and development, human basic necessities of life all cannot do without the land. As land resources are limited, and land resources are non-renewable, but with the rapid increase in population, socioeconomic growth, the destruction of land resources, the conflict between human demand and the limited land resources escalates and the irrational land use has threatened the sustainable development of human society and economy. Therefore, scientific evaluation of the land resources carrying capacity and rational development, utilization and protection of land, integrative distribution of land resources, make sustainable development of land resources possible.

At present, most of scholars at home and abroad concentrate on the exploitation of land resources [START_REF] Yansui | Optimal allocation of regional land use[END_REF], sustainable utilization [START_REF] Hu Yanlan | Method for ecological benefit assessment of rural land consolidation[END_REF] optimal allocation [START_REF] Xiulan | Study on the methods of land use dynamic change research[END_REF], and the ecological environment evaluation [START_REF]Food and Agriculture Organization of the United Nations, FESLM: An International Framework for Evaluating Sustainable Land Management [R[END_REF] and so on, the research on comprehensive evaluation of land resources carrying capacity is very limited. However, comprehensive evaluation is the basis of water-related research aforementioned, Sustainable utilization of land resources can only be achieved by management within its carrying capacity. Land Resources Carrying Capacity (LRCC) is an important indicator for evaluating the regional population, resources and environment and sustainable development [START_REF] Slesser | Carrying Capacity Assessment [R[END_REF], the research has important significance on the development strategy of regional economic development, population growth, land use planning and layout of production. In view of this, land resources not only contains the natural attributes, but also contains the complex social, economic and ecological attributes. Real-coded accelerating genetic algorithm-projection pursuit model was applied to carry on comprehensive evaluation of land resources carrying capacity of Sanjiang Plain, and the comparative study was done under both national and regional scale to form scientific understanding and judgment on land carrying capacity levels. The research results can provide reference for the sustainable use of land resources and the development of society.

optimize the parameters of projection direction in PPC, conversing the highdimensional data to a low dimensional space, that is to say, multiple evaluation indicators are synthesized into an integrated one, then sort by projection value and recognition [START_REF] Qiang | Principle and application of projection pursuit model[END_REF].

PPC Model

PPC model modeling steps are as follows [START_REF] Qiang | Applying PPC model based on RAGA in evaluating soil quantity variation [J][END_REF][START_REF] Jingtian | The statistical method of data processing in nonlinear system[END_REF][START_REF] Xinli | Projection pursuit and its application in hydrology and water resources[END_REF]:

Step 1:

Indicators normalized. The evaluation sample set is   p j n i j i x ,..., 2 , 1 , ,..., 2 , 1 | ) , ( *  
, where

) , ( * j i x
is the jth indicator value of ith sample, n is the number of the samples、p is the number of the indicators. In order to eliminate the variation range of each indicator value dimension and unified the indicator value, the type of normalization, using the following formula for normalized.

The bigger the better indicators (forward):

) ( ) , ( ) , ( max * j x j i x j i x  (1)
The smaller the better indicators (reverse):

) ( ) , ( ) , ( max * j x j i x j i x  (2)
Where

) ( max j x 、 ) ( min j
x are the biggest and smallest value of the jth indicator; And

Where a is the unit vector. Then do the classification according to the one dimensional scatter diagram of  

n i i z ,..., 2 , 1 | ) ( 
. When projecting the indicators, the scattered feature of indicator ) (i z should be as densely as possible for region projection point while as dispersed as possible for the whole projection point. Therefore, the projection indicator function can be expressed as:

z z D S a Q  ) ( (4)
Where z S is the standard deviation of indicator

) (i z , z D is the regional d ensity of ) (i z , that is 1 )) ( ) ( ( 1 2      n z E i z S n i z (5)         n i n j z j i r R u j i r R D 1 1 )) , ( ( )) , ( ( (6)
Where

) (z E is the average value of   n i i z ,..., 2 , 1 | ) ( 
, R is the window radius of the regional density, the usual value is 0.

1 z S ; ) , ( j i r is the distance between samples, | ) ( ) ( | ) , ( j z i z j i r   ; ) (t u is the unit step function, when 0  t , it is 1, when 0  t ,it is 0.
Step 3: Optimizing the projection indicator function. When the sample set of indicators is given, the projection indicator function only varies with projection direction a . The direction of projection reflects the different characteristics of different data structure, the best projection direction projection direction can maximize the feature structure of high dimensional data, so we can estimate the best projection direction through solving the maximum of the projection indicator function, that is:

Objective function:

z z D S a Q   ) ( : max (7) Constraints:    p j j a t s 1 2 1 ) ( : . . (8)
This is a complex nonlinear optimization problem, in which  

p j j a ,..., 2 , 1 | ) ( 
is the decision variables. Traditional optimization method is difficult to handle this problem, so in this paper, accelerated real-coded genetic algorithm is used to address the high dimensional global optimization problems.

Step 4: Comprehensive Evaluation. The best projection direction * a obtained by 3 was put into formula (3) to acquire projection value

) ( * i z
of each sample, i.e. the comprehensive evaluation value of land resources carrying capacity.

Accelerated Real-coded Genetic Algorithm

The genetic algorithm, put forward by Professor Holland of the University of Michigan's America, is a kind of adaptive global optimization probability search algorithm simulating heredity and evolution of biology in the natural environment. It mainly includes selection, crossover and mutation operations. The specific steps are as follows [START_REF] Qiang | Principle and application of projection pursuit model[END_REF][START_REF] Zhou | The principle of genetic algorithm and its application[END_REF]:

Step 1: randomly generate N groups' random variables.

Step 2: calculate the objective function value, sorted by size, from largest to smallest.

Step 3: calculate the evaluation function based on ranking.

Step 4: selection operation, generating new population.

Step 5: crossover the new population generated in step 4.

Step 6: mutate the new population generated in step 5.

Step 7: evolutionary iteration.

Step 8: take the variable interval of outstanding individual generated in step 1 and 2 as the variable initial change interval, then run step 1 step 7, bringing about accelerated process, until the function value of best individual less than the default value or run number reaches the predetermined value. At this time, the best individual in the current population is treated as RAGA results.

Comprehensive evaluation of land resource carrying capacity 4.1 Construction of Evaluation Indicator System

Land resources carrying capacity evaluation indicator system (Fig. 1) was constructed according to the principle of constructing evaluation indicator system, comprehensive analysis on the regional land resources and influence factors, the experience of other scholars' research achievements [START_REF] Shuhua | Design and evaluation on the indicator system of land comprehensive carrying capacity [J][END_REF][START_REF] Ying | Study on the integrated assessment system of land resources carrying capacity for Chongming Island[END_REF][START_REF] Jie | Evaluation of land resources carrying capacity in the process of urbanization[END_REF] and the status of the ecological environment. Evaluation system is divided into the total objective layer, sub objective layer and indicator layer, where the indicator layer contains 29 indicators.

Comprehensive Evaluation of Land Resources Carrying Capacity

In order to be more scientific and reasonable on the evaluation of the land resources carrying capacity in Sanjiang plain, this study selected 10 land samples for carrying capacity evaluation, including the national large-scale, the Heilongjiang Province and Sanjiang plain under medium scale and seven cities in Sanjiang plain under small scale. Value [START_REF]Statistical yearbook of Heilongjiang Province[END_REF]; 2, national evaluation indicator data are in reference [START_REF]The people's Republic of China National Bureau of Statistics[END_REF]; 3, F stands for forward, R stands for reverse.

Land Resources Carrying Capacity Evaluation and Regional Differences

Land resources carrying capacity evaluation on national large-scale shows that in 2012 the national land resource carrying capacity is much lower than that in Heilongjiang province and the Sanjiang plain, the land resources utilization degree is low; Because of population pressure, land resources support capacity for social development and ecological environment protection is low, in the lower level of land resources carrying capacity (grade III), at a high stress level. But the national economic and technological level is high (economic and technical evaluation ranks 1st), hence, we can control the population while increasing funding for land resources for land resources and technology development and management improvement, hence enhance the overall land resources carrying capacity, reducing the social development pressure on land resources.

Land resources carrying capacity evaluation on medium scale shows that because of the rich land resources in Heilongjiang province and Sanjiang plain and the low population density, the evaluation of land resource utilization, the supporting capacity of land resources for social development and ecological environment protection, the comprehensive land resources carrying capacity are higher than the national level, the land resources carrying capacity belongs to middle grade (III) and the use of land resources and social development are at the equilibrium state. Because of the characteristics of geographic location and natural resources, the production activity of Heilongjiang province and Sanjiang plain dominated by agricultural production, industrial structure is irrational, backward in economy. Therefore, in the economic and technical evaluation, the economic and technical level of Heilongjiang province and Sanjiang Plain are lower than the national average level, restricting the advantage role of land resources in this area.

Especially as the commodity grain production base, the Sanjiang plain should strengthen transformation of low yielding fields, put the cultivation techniques of high yield and high yield varieties into practice to comprehensively promote the agricultural production level, improve agricultural land output efficiency, grain production capacity, at the same time, efforts should be made to expand the area of basic farmland protection, improve the role and status of agricultural product base, the formation of high yielding agriculture industry clusters and the industrial chain as soon as possible, hence ensure the use and coordination of economic and technical development of land resources.

Because of the different distribution of natural resource, land resources carrying capacity in Sanjiang plain under small scale varies significantly. The land resources carrying capacity of Jixi City, Hegang City, Shuangyashan city and Jiamusi city are at a high level (grade I), in a stress-free state. The area is rich in land resources, the available land resources quantity is enough to ensure the social and economic development. Future utilization of land resources should be focused on the depth development resources, the rational allocation of land resources, at the same time it should combine with the characteristics of regional resources to adjust the industrial structure, improve the efficiency of resources, to ensure the benign circulation of land resources.

The land resources carrying capacity of Qitaihe City, Muling City and Yilan County in Sanjiang plain belongs to the medium level (grade II),the use of the land resource and social economic development scale is in equilibrium. The area is the area of Sanjiang plain in the highest population density in the area, with the expansion of social development scale, shortage of land resources will become a limiting factor to the social and economic development. So in the future, population density should be controlled while promoting the rational distribution of industrial land, strengthening mining land reclamation, increasing the utilization of barren hills and wasteland, expanding the available land area to reduce the pressure on land resources. The efforts should also be put on the environment protection, the rational development of forest and wetlands, and all the exploitation should be on the premise of ecological environment.

Conclusions

(1) In the full analysis of land resources carrying capacity and regional characteristics, established the regional land resource carrying capacity evaluation indicator system, we set up the regional land resource carrying capacity evaluation indicator system. This system contains one total goal, four sub goals and 29 indicators, evaluating the regional land resource carrying capacity from four perspectives (i.e. land resources utilization, social development, and economy, technological and ecological environment).

(2) The paper constructed comprehensive evaluation model of land resources carrying capacity and acquired the comprehensive evaluation value of land resources carrying capacity under large scale (nation), medium scale (Heilongjiang province, Sanjiang plain), small scale (6 cites and 1 county in Sanjiang plain) in 2012.

(3) After full consideration of the comprehensive evaluation value of land resources carrying capacity, we put land resources carrying capacity into three levels. Among them, the national land resources carrying capacity falls on grade III; Heilongjiang province and Sanjiang plain under medium scale belongs to grade II; Jixi city, Hegang city, Shuangyashan city and Jiamusi city under small scale belong to grade I, while Qitahe city, Muling city and Yilan city belongs to fall on grade II. The result concluded in the paper, coupled with consideration of the characteristics of the regional social and economic development could support the policy making.

Fig. 1 .

 1 Fig.1. Land carrying capacity evaluation indicator system

Table 2 .

 2 table of land resources in different scale statistics 2012 carrying capacity evaluation indicator. The comprehensive evaluation indicator value of land resources carrying different scale in 2012 is shown in Table 1. Matlab is used to program RAGA-PPC model to evaluate the total target and four sub goals of land resource carrying capacity under different scale. The total goal, the best projection value (i.e. Z*) of sub goal and its ranking are shown in Table 2. On the basis of accumulation and dispersion characteristics of land resource carrying capacity value, coupled with the comprehensive evaluation of the land resources carrying capacity in Heilongjiang Province, land resources carrying capacity evaluation level can be divided into three grade: grade I , land resources carrying capacity is weak, in a state of stress. Comprehensive evaluation result of land resource carrying capacity under different scale in 2012 is shown in Table 2. Comprehensive evaluation result of land resource carrying capacity under different scale in 2012

	Z	*  A	30 . 3	, land

*  A Z

Table 1 .

 1 The comprehensive evaluation indicator value of land resources carrying different scale in 2012 the data on cities (counties) in Sanjiang plain and Heilongjiang province evaluation indicators are in reference

	Indicator	Unit	Type	Small scale Jixi Hegang Shuang-yashan Jiamusi Qitahe Muling Yilan Sanjiang Plain Heilongjiang province Medium scale	Large scale China
	C 1	%	F 84.28 90.78	92.86	87.22 96.97 98.44 95.41	89.81	87.1	72.21
	C 2	%	F	79.8 86.52	88.59	82.78 91.41 90.06 90.53	85.11	83.8	68.43
	C 3	%	F 43.74 37.61	48.05	62.91 45.26 23.82 52.89	48.74	26.16	12.68
	C 4	%	F	4.36 4.09	4.16	4.29	5.53	2.01 4.65	4.2	3.3	3.44
	C 5 C 6	% people/km 2	R 0.13 0.33 R 83 74	0.06 68	0.1 74	0.09 149	0.07 0.14 44 88	0.13 78	0.09 85	0.24 141
	C 7	%	F 63.91 80.74	64.3	50.08 58.87 49.48 42.26	60.12	56.9	52.57
	C 8 C 9 C 10 C 11 C 12 C 13 C 14	% hm 2 / capita m 2 / capita m 2 / capita m 2 / capita Yuan/km 2 10 4 Yuan/ capita F R 33.82 F 0.53 0.51 34 F 303 294 F 173 206 F 15.4 26 F 259 244 3.1 3.3	29.8 0.7 311 216 15.3 256 3.7	35.1 38.21 39.61 35.2 0.85 0.3 0.55 0.6 287 222 270 286 287 110 168 197 16.2 26.1 8.2 28.2 206 483 221 260 2.8 3.2 5.1 3	33.35 0.63 288 196 18.5 251 3.5	36.9 0.31 255 195 19.3 303 3.6	36.2 0.09 156 118 12.3 541 3.8
	C 15 C 16 C 17 C 18 C 19	% 10 4 Yuan /km 2 10 4 Yuan /km 2 10 4 Yuan /hm 2 kg/hm 2	F F 92.03 83.15 13.6 13.5 F 4249 4217 F 1.11 0.59 F 3305 2061	13.5 94.65 4148 0.81 2648	13.8 75.04 54.08 43.43 84.95 8.3 17 12.3 3348 7840 9040 3937 1.14 0.84 2.12 1.8 4123 3131 2679 7445	13.2 80.68 4344 1.06 3522	10.1 55.74 7760 1.96 4867	9.7 79.73 14108 3.86 4843
	C 20	Yuan/ capita	F	5.88 10.42	8.56	3.58	2.77	4.53 5.68	5.39	2.65	2.03
	C 21	Yuan/ capita	F 12.12 11.28	18.39	11.55 11.65 26.86 8.24	13.16	16.13	10.12
	C 22 C 23	Yuan/ capita kw/hm 2	F F	5.88 4.58 1.93 1.49	5.18 1.47	6.68 1.65	5.02 2	4.29 5.47 1.48 2.53	5.63 1.7	7.3 3.85	8.36 8.43
	C 24	%	F	56.9 82.2	73.9	50.3	61.8	60.4 80.2	65.8	75.3	51
	C 25	%	F 12.77 19.26	10.44	11.95 39.01 3.57 11.69	13.81	16.31	24.67
	C 26 C 27	% kg/hm 2	F R 46.6 76.1 32.6 45.9	37.7 58.7	16.6 99.2	43.8 48.8	65.7 35.1 58.3 90.8	33.4 74.8	42.6 203	20.4 479.7
	C 28	%	F 11.69 6.8	5.97	10.02 2.55	0.88 4.41	7.89	1.91	1.72
	C 29	%	F 15.59 8.89	7.07	12.69 2.94	1.48 4.45	10.05	9.53	4.01
	Note:1,

Study AreaSanjiang Plain is part of the plains of the Northeast, is China's largest freshwater marsh distribution area, including Ken River, Muling River plains, Xingkai lake plains and the triangular confluence region of Heilongjiang, Songhua and Ussuri. The Sanjiang plain is an important commodity grain production base, and also the key to the realization of Heilongjiang Province billion kilogram of grain production engineering and food safety in China. Over the years, excessive pursuit of economic benefits results in the increasing of plant area in Sanjiang plain[START_REF] Qiuxiang | Spatial matching patterns of land and water resources in Sanjiang plain [J][END_REF], single plant structure. Coupled with the absence of appropriate government, Soil quality degraded severely, resulting in the increase of the proportion of low yielding fields. The utilization of land resources in Sanjiang Plain is regarded to be unreasonable, a series of problems has appeared, such as the low utilization efficiency and the decreasing of the carrying capacity, which has imposed a serious threat to the sustainable use of resources in Sanjiang Plain and sustainable socio-economic development. Therefore, in the premise of ensuring national food security and regional ecological security, scientific evaluation of land resources carrying capacity, the rational allocation of land resources, improving the land resource sustainable support ability in industrial and agricultural production and the development of social economy, seems to be basis of the realization of Sanjiang plain soil resources sustainable utilization.

This study will combine real coding based Accelerating Genetic Algorithm (RAGA) and Projection Pursuit Classification Model (PPC) organically. RAGA is used to

Acknowledgment

The study was supported by the National Natural Science Foundation of China (51179032, 51279031, 51209038), Cultivation Plan of Excellent Talents in New Century of Ministry of Education of Heilongjiang (1155-NCET-004), Ministry of Water Resources' Special Funds for Scientific Research on Public Causes (201301096-0204), Scientific Research Fund of Heilongjiang Provincial Education Department (12531009), Doctor Research Fund of Northeast Agricultural University (2012RCB58) and Heilongjiang Postdoctoral Fund (LBH-Z13049).


Xiaofei An

Zhijun Meng

Guangwei Wu

Jianhua Guo

Software Design of Distribution Map Generation for Soil Parameters Based on VC++

Keywords:

Software system on distribution map generation of soil parameters was developed based on VC++. With the help of MapX control, soil spatial data could be edited, analyzed, generated distribution maps and so on. Soil parameter scatter graph was displayed and classified by equal count method, equal range method, natural break method and standard deviation method. In the software system soil parameter statistics data was also calculated and a new perimeter and area algorithm was developed. After comparison of different interpolation methods, Inverse Distance Weighted (IDW) method was selected to generate the soil parameter contour map and soil parameter iso-surface. At the same time the system gave the legend color and drawing method. These maps could also be saved as vector graphics. Finally a distribution map generation of soil parameters for soil moisture and total nitrogen was conducted by the soil samples from Xiangtang apple orchard in Beijing suburb, including scatter graph, contour map and iso-surface. The errors of perimeter and area were less than 3.72% and 3.01%, respectively. The result shows that the system has perfect function and stable performance.

Introduction

Soil is the basis of agricultural production, human has continued to research and analysis soil in the recent years. According to the theory of soil science, soil has two main functions: the first function is used to store water and supply crops; the second is to maintain the healthy growth of crop root. In order to describe the two functions, human use several parameters such as soil moisture, soil total nitrogen, soil fertility and many other parameters [START_REF] Minzan | Spectral analysis and application[END_REF][START_REF] Lihua | Real-time Sensing of Soil Parameters Based on Spectroscopy[END_REF][START_REF] Jianhua | Research advancement and status on crop nitrogen nutrition diagnosis[END_REF][START_REF] Wangshou | Comprehensive assessment methodology of soil quality under different land use conditions[END_REF] . At the same time, domestic and foreign scholars have obtained a lot of results on the rapid detection of soil parameters [START_REF] Xiwen | Research progress in farming information acquisition technique for precision agr iculture[END_REF][START_REF] Dongjian | Research progress of real-time measurement of soil attributes for precision agriculture[END_REF][START_REF] Mouazen | On-line measurement of some selected soil properties using a VIS-NIR sensor[END_REF][START_REF] Christy | Real-time measurement of soil attributes using on-the-go near infrared reflectance spectroscopy[END_REF][START_REF] Adamchuk | On-the-go soil sensor for precision agriculture[END_REF][START_REF] Minzan | Development of a portable SOM detector based on NIR diffuse reflection[END_REF] . As a kind of new management ideas and modern agriculture, precision agriculture requires precise scientific fertilization and nutrient management [START_REF] Maohua | Precision agriculture" with the practice of agricultural science and technology innovation[END_REF][START_REF] Wu | Research on the variable rate fertilization and economic analysis[END_REF] . In order to achieve the crop production in the process of scientific regulation measures of investment and management decisions, it must implement the precise fertilization decision-making. As a result it is particularly important to generating high precised soil parameter distribution map. Currently, combinations of GIS tool software and visualization technology integrated have become the mainstream of the GIS application [START_REF] Hailiang | Secondary development of GIS application software functions based on MapX controls[END_REF][START_REF] Yong | Information processing system for precision agriculture based on GPS and GIS[END_REF] . This article discussed the development of soil parameter distribution map system based on VC ++. With the help of MapX control [START_REF] Yong | Information processing system for precision agriculture based on GPS and GIS[END_REF][START_REF] Xure | MapX technology based on viscal C++[END_REF][START_REF] Huiwu | The development of GIS Application software on the basis of MapX control in VC++[END_REF] , it could display soil data, scatter plot and generate soil parameter distribution map.

2 The structure of soil parameter distribution map system

Experimental Samples

According to the system requirement, soil parameter distribution generation system mainly included the following modules: basic operation module, data import module, data statistics module and parameter map generation module. Figure 1 is the overall hierarchy diagram of soil parameter distribution generation system. In VC + + 6.0 environment, the system has realized the import of soil parameter information, statistics and mapping, and other functions with MapX control. Soil parameter information included soil moisture value, soil total nitrogen and soil GPS coordinates. Data binding is the process of soil parameter data from the data source into MapX [15]. In the MapX control, every map was corresponded to the multiple layers (the layers), and each layer (layer) have datasets, which contains a dataset object. Datasets with some properties and methods were mainly used to add and remove the dataset object in the collection.

Soil parameter data mainly included GPS location information of soil and soil parameters attribute values. Then the binding layer type was set to miBindLayerTypeXY, which could realize to WGS54 coordinates X and Y coordinates of points. Figure 2 shows the imported the scatter plot after binding the soil moisture data.

基本操作 数据导入 数据统计 分布图生成 帮助

Fig.2 Data binding result

Data classification and statistics

After the soil data was binded to the soil layer, it could be classified by four methods, such as the counting method, the same range method, the nature devision method and the standard devision method. It was said that soil and crop information always has a certain statistical regularity, although sometimes it does not get the distribution function of random variables. In this system the soil parameters could be calculate by maximum, minimum, the statistics of the mean and standard deviation. The mean value and mean square error of the calculation formula were shown in ( 1) and (2).

     1 0 1 n k k v n v (1)       1 0 2 ) v ( 1 n k k v n  (2)
Where:  was the data standard deviation value, v was the average of the data, n was the number of data points, k was the data point number, k v was the first k data of soil parameters.

3.3Perimeter and area calculation

This system designed an algorithm to calculate perimeter and area with the direct use of GPS handheld device data. And projection method was directly to gaussian zoning projection. The center was east longitude 116 °, projection of 40°.

MapX control provided map.distance method. This method could be used to calculate the linear distance between two points. Along the land border manually choose some key points, and then map.distance method was used repeatedly in sequence. Finally the perimeter was obtained accumulation. The system adopts the offset compensation method to calculate polygon area. Figure 4

Soil parameter distribution map generation

With the steps of data binding, scatter plot, interpolation, contour map, soil parameter distribution map could be generated. By comparing the different types of spatial data interpolation methods, inverse distance weighting spatial data interpolation algorithm was finally chosen. After the original data interpolation was completed, the contour map was generated by the contour generating function provided by Contour.OCX control. Both the corresponding smoothness and Contour interval could be set. And finally the map could be saved format of SHP. Figure 5 was the distribution of soil parameters generated dialog.

4.1Efficiency analysis

To test the soil parameter distribution generation system efficiency, different sets of import test data was used. Choose the five sets of data, each group data included the number, soil parameter value latitude and longitude location information. Figure 6 is data import speed test results. It could be seen from the figure 6, when import external data, this system increased with the amount of data. The consumption of time required to import the data were linearly increasing trend. When the data points was 2000 sets (the amount of data with 8000), the system took 1344 ms. From the perspective of the function of angle and speed of system, functions could be fully realized, moderate speed in data quantity could meet the needs of practical application.

Perimeter and area validation

In order to test this system the perimeter area of the accuracy of the algorithm, experiments were carried out in Beijing changping apple orchards with GPS handheld instrument G738CM and the system. After sample the apple orchard with the device, the perimeter and area were calculated by the system. Table 1 It could be seen from table 1 and table 2, with the increase of the number of sampling points from 6 to 28, apple orchards perimeter error reduced from 3.72% to 1.42%, the average was 2.46%; Area error reduced from 3.01% to 2.04% , an average of 2.14%. Key points directly affect the precision of the result. The more the key point was, the higher precision was. The error of the system was less than 3.72% and 3.01%, respectively. It could satisfy practice.

Soil parameter distribution map validation

In order to verify the system operation, soil moisture and total nitrogen distribution maps of changing apple orchard were generated respectively, including the scatter plot, contour map and level surface. Inverse distance weighting spatial data interpolation method was used on the original data interpolation. The results showed that the system had stable performance. Soil parameter distribution map system was developed based on VC++. The system could realized the scatter plot, contour map, contour surface generation and vector graph save function. Through the analysis of the system efficiency and perimeter area computation verification test it showed that when the data points was 2000, the timeconsuming needs 1344 ms, calculation error within 3.72% and 3.01% respectively; The results showed that the system was stable and reliable enough. It provided a new means to guide the variable rate fertilization.

  Fig.1 Overall structure of system

Figure 3 (Fig. 3

 33 Fig.3 Data classification method

Fig. 4

 4 Fig.4 Area calculation flow chart

Fig. 5

 5 Fig.5 Soil parameter distribution map

Fig. 6

 6 Fig.6 Data import speed result

Figure 7 (

 7 a), (b), (c) and figure 8 (a), (b), (c) were the distributions of soil moisture and soil total nitrogen.

  and table 2 were the perimeter and area calculation error, respectively.

	4	20	245.93	241.06	2.02
	5	28	245.58	242.14	1.42
			Tab.2 area calculation error
	Number	Key point	System area/m 2	System area with	error/%
				GPS device /m 2
	1	6	2160.62	2097.40	3.01
	2	8	2079.50	2035.34	2.17
	3	12	2110.55	2087.28	1.11
	4	20	2231.84	2180.40	2.36
	5	28	2248.96	2203.98	2.04
			Tab.1 Perimeter calculation error
	Number	Key point	System	System perimeter	error/%
			perimeter/m	with GPS device/m
	1	6	247.00	238.14	3.72
	2	8	246.81	238.00	3.70
	3	12	242.90	239.47	1.43

Acknowledgment

This paper was supported by the 863 programmar (2012AA101901) and 948 programmar(2014-G32).