Mingxuan Li

Man Zhang

Haiyan Huan

Gang Liu

A new navigation line extraction method for agriculture implements guidance system

Keywords: Navigation line detection, Machine vision, Agricultural implements navigation

To make up the shortages of the existing vision navigation algorithms such as the noise interference and the processing speed, a new approach for navigation line detection is presented in this article. This new method was based on image boundary detection. The color images were transformed into binary ones after pre-processing. Then, three stages in this method were presented. At first, the image was divided into several sectors and the boundary of the crops was detected. Then a number of points indicating the centers of the rows were determined. At last, the navigation line was estimated by an improved linear regression algorithm. Compared to other navigation line detection algorithms, this new method could find the navigation lines more accurately. It need only about 89ms to process a picture with 640×320. So it can meet the need of agricultural real-time visual navigation when the speed is 1~2m/s.

Introduction

Steering agricultural machinery within rowed crop fields is a tedious task for producers. Automated guidance of the machinery will not only reduce operator fatigue but also increase both the productivity and safety of the operation [START_REF] Marchant | Real time tracking of plant rows using a Hough transform[END_REF] . As a machine vision system can detect a path in relation to crop rows, vision-based guidance can be used to guide machinery travelling between crop rows to perform field operations such as cultivating, chemical spraying, and harvesting [2~4] .

There are a lot of research have been made in the field of crop line identification. At present, the traditional ways of crop detection method is based on the Hough transformation largely [5~7] . Hough transformation algorithm is a shape matching technology put forward by Paul Hough in 1962. Duda and Hart established the Hough transformation linear testing in 1972. Hough transformation has several advantages such as strong robustness and small noise influence [START_REF] Jang | Research on Navigation and Localization of Agricultural Machinery Based on Machine Vision[END_REF] . However, it also has disadvantages such as a high computational complexity [START_REF] Cao | Detection algorithm for crop target multi-lines of the field image based on machine vision[END_REF] . In order to improve the speed of Hough, Kultanen P. put forward the Random Hough Transformation (RHT) [6, 10] algorithm. RHT is a technology of multi-to-one. Compared to the traditional Hough transformation, RHT depend on the idea of randomness, reducing the computational complexity greatly. But for the RHT algorithm, random point selection is so important that it is unstable to detect the navigation line [START_REF] Liu | Methods of Navigation Control and Path Recognition in Autonomous Agricultural Vehicle[END_REF] . In addition, because the navigation line only depends on the two crops lines nearest to the image center, so it still brings some redundant calculation to extract all the crop lines.

In order to solve these questions above, a new method was introduced in this article. Multistage filtering is put forward to solve the detection problem, and result in a lower computational complexity with the similar detection stability than RHT and traditional Hough transformation. Two sections are contained in this new method, the first one is image pre-processing section and the second one is a new baseline extraction algorithm that designed in this research named Boundary Detection Navigation Algorithm.

2

Materials and method

Experimental Materials

The experimental platform, as shown in Fig. 1, included several parts, such as vision sensor, computer, controller, hydraulic system and the agricultural implement. In this platform, the camera was used to real-time acquire the image of farmland; the computer was used to identify and position of crop line from the acquired image, and then calculate the lateral deviation of the farm machinery relative to the target path; Based on the lateral deviation, the controller is used to control the hydraulic system work and realize the automatic navigation of fame machinery.

Fig. 1. The experiment platform

The visual sensor used in this paper is OK AC1310 camera made by Beijing Jiaheng Company. The pixel size is 4.65 μm×4.65μm; the resolution is 1300 × 1024; The frame rate is between 8~24Hz; The shortest exposure time of the electronic shutter is 20μs. The data transmission between the video camera and the computer through a USB2.0 interface [START_REF] Chen | Crop rows detection based on parallel characteristic of crop rows using visual navigation[END_REF] .

The navigation program is developed on the computer of Intel(R) Core (TM2) Quad CPU Q9400 @2.66GHz ,2.96GB in the XP operating system. And the development environment is Microsoft Visual C++ 6.0.

Methods

The testing method of the navigation line was divided two sections, the preprocessing section and the navigation line extraction section. The pre-processing section included graying, binarization and filtering, and the navigation line extraction section included division and detection (including image segmentation and boundary detection), points determination (including false points removing and navigation points estimation), and navigation line regression, as shown in the Fig. 2. Considering the complexity of the farmland environment, color images were chosen as research objects in this article because rich color information was contained in it [START_REF] Liu | Methods of Navigation Control and Path Recognition in Autonomous Agricultural Vehicle[END_REF] . However, color images also brought some problems such as the redundant information, the noise distribution and the complex nature. Thus, the images needed to be preprocessed before navigation line extraction. The image preprocessing contained three segments, namely as graying, binarization, and de-noising.

Usually, the green component G is much greater than red R and blue B components because the main pigment of the crops was chlorophyll. So the 2GRB method was taken to gray images in order to emphasize green component and restrain the rest two components [START_REF] Chen | Crop rows detection based on parallel characteristic of crop rows using visual navigation[END_REF] . The brightness and color will different if the image is acquired under different light condition. The OTSU method was used to transform the gray images into binary images [START_REF] Soggaard | Determination of crop rows by image analysis without segmentation[END_REF].

There were a lot of noise points in the image after segmentation because of the characteristics of farmland environment, such as weeds, shadow. The method of corrosion-median filter-expansion was selected to reduce the noise according to the characteristics of the weeds. Firstly, the disc structure with 5 pixel radius was chosen to corrode the images. Secondly, the image was filtered through the 5×5 pixel size windows. Finally, the 5×1pixel inflation structure elements was adopted to extrude the vertical part of the images. After that, most of the noise was filtered out and the stability of the navigation baseline extraction algorithm would be improved.

The images of bean, corn and winter wheat were collected to verify the effect of the image preprocessing. The results of the bean images were shown in the Fig. 3, it showed that the preprocessing algorithm could effectively reduce the influence of weeds in images and get the boundary of the crops.

The navigation line detection algorithm

Before detecting the row position, a number of points indicating the centre of the rows needed to be determined firstly. In order to obtain these points, the binary image resulting from grayscale image was divided into a series of horizontal strips as shown in Fig. 4, and the next step was to estimate where the rows were intersect each individual strip. In principle, the number of strips should equal to the number of pixel rows in the image. In order to reduce the amount of subsequent computation and enhance the stability of subsequent analysis, it was suggested to let each strip involve more than one pixel row. Thus, the number of strips should be kept at a reasonably low value, and this value was denoted as step . Practical experience showed that every 10 pixels divided into a strip was suitable when the system was working on the experimental field. Each image strip was processed as illustrated in Fig. 4. Mathematically, the process could be described as Fig. 5. was used as a row vector resulting from summing the array V over the index i:

, 1 N i i j j V     1, 2... iM  (1)
The result of this summation is illustrated by the curve in Fig. 7. The position of non-zero value in Fig. 7 presents the position of crop in Fig. 6. So. Crop position can be determined by detect the point that exchange between 0 and 1. The ith non-zero part is represented by a 2D vector like P i =(Pl i ,Pr i ). Pl i was the position of the ith nonzero part's left border in a strip, Pr i is the position of the ith non-zero part's right border in a strip. However, not all of these non-zero points represents the real crop because of the influence of weeds. An important difference between crop and weed is their scale in the image. By make use of this feature, a criterion to distinguish the crops and weeds could be set. This criterion describe as follows:

The ith non-zero part represent the crop when

1 1 ( ) ( ) 2 K i i j j j Pl Pr Pl Pr K     
, (K is the number of non-zero parts in a strip),or it represent the weed.

Step 4: Navigation Points Estimation

After previous step the navigation points ( , )

navx navy ii PP can be estimated by the following formula:

0 1 min (| ( ) |) min (| ( ) |) 2 2 2 2 2 * 2 i r i r navx i navy i

Pl Pl

Pl Pl width width

P step P step i               (2)
In this formula, the 0 min ( ) S is the smallest value in S , the 1 min ( ) S is the second smallest value in S , 1, 2... iN 

, and the navigation point detection result is shown by Fig. 8.

Fig. 8. The detected Navigation Points

Step 5: Navigation Line Regression

Based on these navigation points, the navigation line can be calculated by LSE (least-square estimation) algorithm. However, because of the existing of surplus plant or lack of one or more plant in a row, there are some outliers exist in navigation points, and not all of navigation points should be applied by LSE algorithm. We set another criterion to determine whether a navigation point is a outliers as follows:

The ith navigation point is a outliers when var( )

navx navx i PP   , and
var( ) P is the variance of P . In our experiment the  is 1.5.

Assume the remainder navigation points after removed the outliers are

           
(3) Where:

22 1 1 1 1 1 1 1 1 , , , ( ) ( )
P N N N N            
The navigation line detection result is shown in Fig. 9.

Fig. 9. The detected Navigation Line

Evaluation test and dicussion

This algorithm was named Boundary Detection (BD) algorithm because it partly depended on crop boundary detection. In order to evaluate the speed and accuracy of this BD algorithm, all traditional Hough Transform (HT) algorithm, Random Hough Transform (RHT) algorithm, and BD algorithm were used to analyze the wheat images. The aim of this test was to extract the center line of the image and calculate the parameters of the navigation line. The result was shown in Fig. 10. The line parameter results were shown in Table 1. These lines were represented by the polar coordinate, and the ρ resented the polar diameter, the θ represented the polar angle. According the above analysis, all of these methods could detect the navigation line successfully with acceptable error. From the timeconsuming perspective, HT consumed the longest time and RHT algorithm was the second. BD algorithm needed the shortest time, which was 89ms and less than half of the HT. From the precision perspective, the difference of θ value between BD and HT algorithm was 0.2567°, and 0.43° between BD and RHT algorithm. And the difference of ρ value between BD and HT algorithm was 2.26mm, and 1.88mm between RD and RHT algorithm. The accuracy of all these methods could meet the demand of the agricultural navigation, but the BD algorithm required the least time.

Conclusion

A new navigation line extraction algorithm named BD algorithm was proposed to process the farmland environment images in this article. After images acquisition, firstly the color images were transformed into grayscale images by 2G-R-B method according to the characteristics of the farmland crops. Then the grayscale images were transformed into binary images by the OTSU method, meanwhile, a corrosionmedian-expansion filter algorithm was adopted to filter the noise. Finally, the parameters of the navigation line were calculated by the BD algorithm. Theoretical analysis and experimental test results showed that the designed BD algorithm was more efficient than HT and RHT algorithm, which only needed 89ms to process an image in this article. The accuracy of the detected line was higher and could be used as the navigation baseline.

Fig. 2 .

 2 Fig. 2. The flow chart

Fig. 3 .

 3 The results of image processing

Fig. 4 .Fig. 5 .

 45 Fig. 4. Color image divided into horizontal strips

Fig. 6 .Fig. 7 .

 67 Fig. 6. The first strip of binary image

Fig. 10 .

 10 Fig. 10. The results of navigation line detection

Table 1 .

 1 The results of linear parameter comparison.

	Time

Method

 R

Acknowledgements

This work was supported by the National Natural Science Fund (Grant No. 30900869), National Technology Support Program (Grant No. 2011BAD20B06), and 863 Program(Grant No.2012AA101901).


Jianlin Li

Lei Cong

Yongxia Yang
email: yongxia.yang@qq.com

Dehai Zhu

Min Yin

Definition and Standardization of Data Elements' Attributes in Land and Resources Management

Keywords: Standardization, Data Element, land and resources, attribute

ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

The standardization of data element is an important component of the standardization of information technology. It guarantees the development and running of information system. The success of the system applications is on the basis of it.

The standardization of data element is defined as the process of developing and implementing uniform standards, which includes the data elements' general rules, definition, description, classification, expression and registration.

Its importance can be broadly divided into five areas:  Standardizing and unifying the collection and application of data;

 Building the basis of data sharing for different systems with multiple modes of data;  Achieving data management in the level of metadata and making the separation of software and data, so the development of software is more flexible.  Providing unified standard of data exchanging in business internal areas or between internal and external areas;  Being the basis of building a unified, integrated and efficient data model in business.

Land and resources management is a multi-disciplinary integrated system relating to land, geography, minerals, etc. Business management needs cooperation and mutual penetration in the area, while it actually lacks of common data models and standards. The present situation causes many problems, such as one word with several equivalents of different meanings and one word with multiple equivalents of the same meaning. The former can make users understand the word with ambiguity like Digital Gap; the later can cause data redundancy and update data inconsistently, and further, it goes ill with the integration of information system, which is called Information Island.

These problems add the confusion of the system, and limit data exchanging and sharing. In summary, the standardization of data element in land and resources management is a vital foundation in this field, especially for the goals of "A Map" and "Golden Land Project". This study is on the standardization of data element in Beijing Municipal Land and Resources Management, especially cadastral management and land use management. The above standards are the basis of the study, whereas, because of language, expression, area, data and application, there are still some problems that need further research. In this paper, we record and summarize the issues and our solutions that focus on the standardization of data element's attributes.

ISO (the

Definition of Data Element' Attributes

Data element is described by a set of attributes. The basic attributes of data element can be organized into five categories, including a total of 22 attributes as following [3] . The above attributes are just recommended in the standards. In the actual study, they can be added, deleted and modified with many factors. First, we can identify priorities on attributes considering of the application of data elements. For example, if the standardization of data element is for the internal needs of an organization, the attributes of Administrative should be completed and detailed for future use in business management; while the standardization of data element is for data sharing among different organizations, the attributes of Administrative may have no much sense. Second, because of the difference of languages and areas, we can add some specialized attributes adapt to the current environment, such as Chinese name, English name, Chinese spelling. In addition, we can also merge and resolve some attributes in order to make them clearer and more concise.


Considering of the actual need in business areas of Beijing Municipal Land and

Resources Management, we determine the following factors to choose the attributes with the aim of less attributes but also a clear expression.

1) Requires attributes, such as name, identifier, definition, data type, etc.

2) Attributes for management, such as version, registration authority, submitting organization.

3) Attributes for inquiring and reading by user, such as keywords, synonymous name.

By the above factors, we define the attributes which are used to describe data elements in land and resources management (see table 1).

Table1. Attributes of Data Element in

Standardization of Data Element's Attributes

This section shows the key points of the standardization of data element's attributes in our study fields, which are cadastral Management and land use management. We study and develop many proprietary and practical rules to solve the problems caused by language, expression, area, data and application. This study complements and improves the processes and methods for the standardization of data element's attributes in specific field, which are abstract and conceptual in the standards.

Chinese Name

Chinese name is the name of data element with single or multiple Chinese words.

It's generally composed of object class term, property term, representation term and sometimes also qualifier term, which are related to the data model of data element. In order to make data element's concept, classification and representation more clearly, we add the first three terms in the attributes of data element, and the methods to identify them are showed in the ensuring chapters. In this part, the development of naming rules is the key, even in our whole research. A standardized name is the basis to solve "Digital Gap" and "Information Island" in the system.

Besides the basic naming rules in the standards as following [5] [7] :

Semantic rules:

 One and only one object class term shall be present.

 One and only one property term shall be present.

 One and only one representation term shall be present.

 Qualifiers will be added as needed to make the name unique within a specified context. The order of the qualifier terms is not significant. The property term is "姓名" ,and the representation term is "名称". They make the meaning repeated, so the representation term is deleted in the name.

 Qualifier terms shall be prefixed to the part qualified. The order of qualifiers shall not be used to differentiate names.

EXAMPLE

"变更前宗地地籍号" and "宗地变更前地籍号" both mean parcel number before changing. The position of Qualifier term "变更前(Before Changing)" have no effect.

Uniqueness rule:

All names in each language shall be unique within one context.

We also develop the additional lexical rules to adapt to the particularity of this field, considering of the use of Chinese words.

Lexical rules:

 These are all of rules we have studied and developed for naming data elements.

As there aren't unified English name for data elements in this field, we don't consider English name as one attribute of data element.

Chinese Spelling

Chinese spelling is a specific attribute for Chinese data element, which is composed of the spellings of Chinese characters in the name of data element. All the letters should be lowercase, and connect with a hyphen"-".

There are two points needed to pay attention to: NOTE 1:

"u" and "v" should be distinguished properly and used in the right place, such as "吕(lv)" and "区(qu)".

NOTE 2:

If the Chinese name of data element with Arabic numerals, the Chinese spelling should retain the Arabic numerals.

For example,

线 状 地 物 权 属 单 位
1(xian-zhuang-di-wu-quan-shu-dan-wei-ming-cheng-1 ).

Identifier

Each data element should have an identifier, which is defined by the registration authority. This process includes the study of classification of data elements. There we suggest that the authority just provides the idea of classification instead of changeless classification number, because the standardization of data elements is an increasing long-term work.

Since there are so many horizontal and vertical business areas in land and resources management, we use the method of classification combining surface and line. First, we divide the field of land and resources management to three areas of land, resource and geological environment by the means of surface classification. Second, we use the combinative method to get business class and business sub-class. After that, the classification can be refined by the object class that data elements belong to. Fig. 1 shows the classification result of land and resources management.

Then we can define identifier of data element based on the above result. Fig. 2 shows the frame of identifier. 2) Business Code is the code of business types. It is represented by two numbers.

3) First-level Classification Code is about the object class that data elements belong to. Third-level Classification Code: 001

Synonymous Name and Context

Synonymous name refers to the different names in the different application environment for one data element. A data element may have none, one or more synonymous name, which are separated by semicolon. The collection of synonymous name can further solve "Digital Gap" and "Information Island" in the system, considering of the traditional habits in the business management.

As synonymous name's corresponding, context records the application environment and process generating and using the data element. A data element may have none, one or more context, which are separated by semicolon, that's decided by the quantity of synonymous names.

Definition

Definition is the statements used to express the essential characteristics of data element. It can distinguish data element from all other ones. We should make them brief and accurate. Besides, definitions should be unique in a field and use singular form. Reading various terminology standards for reference is the best way to make definition.

Object Class Term

Object class term can show the objects and concepts that data elements belong to, In order to prevent object class term confused with qualifier term, the entities and relations must be minimal. The object class that data element directly belongs to, which isn't able to be subdivided, generates the object class term of data element, while the object class at the higher level may be the qualifier term. Sometimes the object class is too small to present in the ERD, so it uses the object class at the higher level instead of the direct one.

and

Property Term

A property is a characteristic common to all members of an object class. Each property has a name. Property term represents the characteristics to distinguish the object classes of data elements significantly. It generally corresponds to the attributes in the ERD, such as address, name, price, etc.

Representation Term

Representation term describes the form of the set of valid values of an administered item which includes representation. Often, the representation term may be redundant with part of the property term. When this occurs, one term or part of one term may be eliminated in a structured name. This can be established as a naming rule.

Representation term is helpful to achieve the functions of data element naming, analysis, analog and inquiry. The international used representation terms include date, code, name, sum and so on. They can also be added by the demand in the study field.

Relationship

Relationship describes the connection between the present data element and other ones. We can know the deep-rooted network of relationship in this business area by the attribute. It includes the following basic relationships: derive-from, compose-of, replace-of, link-with [START_REF]GB/T 19488 Date Elements for e-government[END_REF]. we rule that the data format of spatial data is represented by "O", and the data format of binary data is represented by "B".

EXAMPLE

Value Domain and Unit of Measurement

Valued domain shows the allowable value set of data element. It's given by four ways as following [START_REF]GB/T 19488 Date Elements for e-government[END_REF] :

1) By name, pointing out the name of valued domain directly.

2) By reference.

3) By list, representing every possible value with its instance or meaning.

4) By rules indirectly.

Unit of measurement should refer to the related international standards. It needs to be represented in Chinese to avoid users' not understanding the abbreviations.

Status

Status shows the current state of data element for registration during its whole lifetime. We define four types of status for data element in land and resources management.

1) Draft: the relevant departments and units can give widely comments and suggestions during this stage.

2) Trial: the result can be used for a try after a certain technical procedures.

The feedbacks in this stage are some practical advices.

3) Implementation: the content of data element has become the standards to be complied with.

4) Abolished: the content of data element would be deleted from the standards.

Discussion and Conclusion

This study is based on the theories and concepts provided by the international standards about the standardization of data element; in addition, we consider of the particularity of data in land and resources management from the aspects of language, expression, area, data and application. The proper attributes of data element in this field are defined and the specific rules to describe them are developed in our study.

They can direct the other related working in land and resources management, and give the ideas to carry out the study on the standardization of data element in the other fields.

Although the study has been basically completed, there are some uncertain issues needed for further study.

1) Chinese name is too long. That's because of the semantic rules. The rules make the name standardized, but also different from the previous used name.

Can we add an attribute called "short title" to record the previous used name?

2) Synonymous name and context are two corresponding attributes. Context

shows the synonymous name's application environment. How to show the application environment of Chinese name?

3) We define the data type and data format of spatial data. But they can be only used for small range. There should be an international standard to define them.

By this study, we have achieved more than 400 data elements' standardization in the areas of cadastral management and land use management. We complete and improve the whole process and rules during the case study. Besides, our research results also include a work guideline to direct the related work in the future. Further, we are making an effort to develop it up to a local standard in Beijing. For the unified management of the standardized data element in the future, a management information system is needed. It ensures the interoperability of data elements in the different departments and business, and achieves data sharing successfully. So the design, development and maintenance of the MIS will be the important work.

This study is only about two areas of cadastral management and land use management. There may be still many problems in the case study of other fields. So the results of this study needs to be improved constantly.

  International Organization for Standardization) and IEC (techology-Specification and standardization of data elements" in 1999, and updated it to "Information technology--Metadata registries (MDR)" in 2004. It's an international standard for defining, describing, standardizing and registering data element or metadata for an organization. For the standard's importance, China National Institute of Standardization translated it into "GB/T 18391-Information techology-Specification and standardization of data elements" in 2002 and "GB/T 18391-Information technology--Metadata registries (MDR)" in 2009 for Chinese organizations' use. In addition, CNIS also released "GB/T 19488 Date Elements for e-government" for better understanding and more practical, which refines and embodies the contents of principle and theories in GB/T 18391.

Fig. 1 .Fig. 2 .

 12 Fig.1. The classification result of land and resources management

  represent an activity or object in a context. The extraction of object class term is a problem in our study, because it usually confuses with qualifier term. To solve this problem, we introduce the thought of modeling, such as Entity Relationship Diagram (ERD). an object class is a set of ideas, abstractions or things in the real world that are identified with explicit boundaries and meaning, and whose properties and behavior follow the same rules. It corresponds to the entity or relation in the data conceptual model. Through rebuilding the ERD, We can sort out all the entities and relations in the whole business areas. The results are the object class and their attributes are data elements.

Parcel compose-of Boundary lines 3 . 10 Data

 310 Type and Date Format Data type is the type of data element value. The basic data type includes character, number, Date, Date Time, Boolean and binary-type. Because of the particularity of data in land and resources management, we should define the additional data types for spatial data. Based on the standards of OPENGIS, we add seven data types of spatial data, which are point, linestring, polygon, multipoint, multilinestring, multipolgon and geometrycollection. Data format orders the format requirement for data elements in a business perspective. It combines three basic attributes which are data type of data element values, maximum size of data element values and minimum size of data element values. There are the exact basic requirements in the existing standards. Besides them,

  If the object class of data element is human, its name should use "姓名" to represent. While the object class of data element is thing or organization, its name should use "名称" to represent.Although they are the abbreviations and acronyms of English words, almost everyone knows and consistently uses them in this field. On the contrary, if they are translated to Chinese, someone may have doubts about them and the name of data element would be too long. Arabic numerals can be used in the name to distinguish data elements.

	than one word can correspond to one meaning but also have a little difference. So
	rules are much needed to prevent the confusion.
	 International abbreviations and acronyms are allowed to use in the name.
	EXAMPLE
	DEM	Digital Elevation Model
	DTM	Digital Terrain Model
	EXAMPLE
	宗地法人代表姓名 Name of Parcel Legal Representative
	宗地权属单位名称 Name of Parcel Ownership Unit
	Comparing to a word of name in English, it seems complex in Chinese. More

EXAMPLE

线状地物权属单位 1 Name of Linear Feature Ownership Unit 线状地物权属单位 2 Name of Linear Feature other Ownership Unit This is a means to make the name of data element simple and brief, but it's on the premise of a clear definition.

Acknowledgements. This research was funded by China Agricultural University

Young Teacher Innovative Project "study on the standardization of data element in Beijing Municipal Land and Resources Management". Young Teacher directed and contributed to the whole study. The authors would like to thank the anonymous reviewers for their constructive comments and suggestions.

Code Business Code First-level Classification Code Second-level Classification Code Third-level Classification Code Land and Resources Management


Sai Tang

Jianyu Yang
email: ycjyyang@cau.edu.cn

Chao Zhang

Dehai Zhu

Wenju Yun

Design of Monitoring Network for Cultivated Land Quality in County Area Based on Kriging Estimation Variance

Keywords: cultivated land quality, spatial sampling design, monitoring network, Kriging, estimation variance 1

Through the monitoring network for cultivated land quality in county area, distribution and changing trend of the quality should be reflected. Besides, the quality of non-sampled locations should also be estimated with the data of sampling points. Due to the correlation among spatial samples, traditional methods such as simple random sampling, stratified sampling and systematic sampling are inefficient to accomplish the task above. Thus, spatial sampling method based on Kriging estimation variance is presented in this paper. The method considers the spatial structure characteristics of cultivated land and is configurable to restrict the expense or estimation accuracy, making it more flexible to determine the number of samples. What's more, tessellation polygon is used to densify the network and this helps Kriging interpolation to achieve higher efficiency than traditional methods. Furthermore, inference for the population mean of natural quality indices stays in a very precise level at the same time.

Introduction

China is an agricultural country with a large population but not enough cultivated land. Until 2011, the cultivated land per capita in China was 1.38 mu (0.09ha), only 40% of the world average. With the rapid development of economy, industrialization and urbanization, this situation becomes even worse [START_REF] Shi Shuqin | Methodology for Impact Assessment of Regional Cultivated Land Resources Change Based on 3S Technology[END_REF]. Therefore, the rate of change of cultivated land is concerns for Chinese governments. At present, methods for monitoring the quantity of cultivated land have become reliable; some researchers analyzed and interpreted remote sensing images to keep track of the quantity condition [START_REF] Jianhong | Accuracy and Efficiency of Different Spatial Sampling Schemes for Cropland Change Monitoring[END_REF][START_REF] Jinshui | Analysis of Influence Factors about Space Sampling Efficiency of Winter Wheat Planting Area[END_REF][START_REF] Hu Tangao | Researches on Remote Sensing Area Measurement Based on Different Sampling Methods[END_REF][START_REF] Xianfeng | Paddyrice Area Estimation Using a Stratified Sampling Method with Remote Sensing in China[END_REF][START_REF] Wang | Spatial Sampling Design for Monitoring the Area of Cultivated Land[END_REF]. However, unreasonable land exploitation in China leaves a negative influence on both quantity and quality of cultivated land. When tracking quality condition, the method of remote sensing is not applicable any more so that on-the-spot sampling and investigation become necessary.

As the dynamic monitoring of cultivated land quality is the next task of Ministry of Land and Resources [START_REF] Yupeng | Research on Standard-plot Based Monitoring and Early-warning of Arable Land Quality[END_REF], we have sampled seven hundred counties in China and in each county, a network would be set up to monitor the change of cultivated land quality. With the basis of current Land Resource Survey and Cultivated Land Evaluation, the paper centers on how to design the monitoring network in county area.

Study Area and Data

The study area Daxing district, located in the south of Beijing (longitude 116°13′-116°43′E, latitude 39°26′-39°51′ N), is a river alluvium plain with a sub-humid warm temperate continental monsoon climate. The terrain slopes gradually from east to west and the span of altitudes is generally between 14 and 52 meters. The average annual temperature and precipitation are 11.6℃ and 556mm.

The data used in our study was the layer of cultivated land quality of Daxing district from last investigation and the corresponding scale was 1: 10000. It included 5182 map spots, with the total area of 416.97km 2 . Daxing district belongs to the region of Huanghuanghai and the sub-region of piedmont plain of Yanshan and Taihang Mountains [START_REF]Ministry of Land and Resources: Regulations for Classification on Agricultural Land (TD/T1004-2003)[END_REF]. When designing the network, only natural quality was concerned and natural quality index was used as the attribute variable. The reasons not to take utilization and economy quality into

Methods

The design of monitoring network is to sample from 5182 blocks of cultivated land and place monitoring points on each sampled block, then infer the quality of non-sampled blocks based on investigation of monitoring points. Currently, most designs are originated from classical sampling models, which mainly include simple random sampling, stratified sampling, systematic sampling and cluster sampling [START_REF] Kish | Survey Sampling[END_REF].These methods ignore the spatial correlation [START_REF] Jinfeng | Spatial Analysis[END_REF] so that they are not suitable for our study because cultivated land belongs to spatially distributed object and there are correlations among samples. Therefore, sampling method based on Kriging estimation variance is presented in this paper.

Kriging method was first raised by Matheron [START_REF] Matheron | The Intrinsic Random Functions and Their Applications[END_REF][START_REF] Matheron | Kriging or Polynomial Interpolation Procedures --a Contribution to Polemics in Mathematical Geology[END_REF], originally applied to researches of geology mineral resources. Oliver used it in the interpolation of GIS [START_REF] Oliver | Kriging: a Method of Interpolation for Geographical Information Systems[END_REF] and later on the method was widely spread in many fields. It is built upon the theoretical analysis of semi-variance function, considering the distance between non-sampled and sampled points. Kriging method belongs to optimal unbiased estimation so its prediction is better than that of traditional interpolation methods. By interpolating the data of monitoring points, we can get the relationship among estimation variance, sampling distribution, spatial correlation and number of samples [START_REF] Jinfeng | Spatial Sampling Design for Monitoring Drought, Flood and Earthquake Disaster in China[END_REF].

Examination of the data distribution

Measures of monitoring network performance

The performance of a monitoring network can be measured using different indices. The following lists three common measures [START_REF] Wang | Spatial Sampling Design for Monitoring the Area of Cultivated Land[END_REF]:

Measure 1: MSE (mean square error) of the mean estimation

    2 n 1 i 2 n A 1 n 1            A i i i dx x f x f E  (1)
Where E denotes the mathematical expectation; n is the total number of samples;

  i x f
is the real value of point x i in the region A; and term

   A i i dx x f A 1
is the true mean of A. In practice, this mean is estimated by the arithmetic mean of n sampling points and the corresponding formula is

    n 1 i n 1 i x f
. This measure is quite useful for estimating population mean [START_REF] Griffith | Heterogeneity of Attribute Sampling Error in Spatial Data Sets[END_REF] or total [START_REF]Elements of Sampling Theory[END_REF].

Measure 2: MSE of discrepancy                A f f E f MSE i i i 2 x - x x ˆ (2)
Where   i f x denotes the real value at point

i x and   i f x ˆ is the estimated value of   i f
x , which is calculated from interpolation of monitoring points. This measure is suitable to estimate other points by spatial interpolation on existing sampling points [START_REF] Griffith | Statistical Analysis of Spatial Data in the Presence of Missing Observations: a Methodological Guide and an Application to Urban Census Data[END_REF].

Measure 3: Maximum site error             A i f f i i i , x - x max (3)
The applicable condition of this measure is similar to that of measure 2.

The quality of non-sampled cultivated land should be predicted by interpolation in our study. As administrators only concern the global accuracy so that measure 2 is appropriate, which places sampling points in suitable positions and has the most accurate estimation through Kriging interpolation [START_REF] Simbahan | Sampling Optimization Based on Secondary Information and Its Utilization in Soil Carbon Mapping[END_REF]. Because natural quality was the prior knowledge, formula of measure 2 was redefined like this:

         n 1 i 2 x - x n 1 i i f f MSE (4)
Where n means the total of cultivated land map spots, which was 5182. By the definition, the higher the estimation accuracy of interpolation is, the lower the MSE is.

3.3

Design of monitoring network

The steps to set up the monitoring network are as follows: a) Lay out the basic network according to regular grids; the choice of grid size depends on both estimation variance and number of sampling points. b) Adjust the centers of grids which do not fall in the area of cultivated land to form the global monitoring network. c) Densify the global network by adding sampling points in tessellation polygons whose MSEs are relatively high. The workflow is displayed in figure 3.

Basic monitoring network.

In geographical space, using regular grids to place samples is the common way to quantitatively describe spatial variation and a square grid is the most convenient in practical applications [START_REF] Renzhong | Spatial Analysis[END_REF]. The reasons to introduce regular grids are listed in the following: a) Sampling points are distributed uniformly so that the entire district could be represented. b) When applying Kriging interpolation, if numbers of point pairs of different lag classes vary greatly, the estimation efficiency of semi-variation function will fall [START_REF] Chengsheng | A Review of Geo-Spatial Sampling Theory[END_REF]. c) Generally, spatial correlation declines when the distance between monitored objects increases so that the grid sampling is more effective than random sampling in decreasing estimation variation [START_REF] Jinfeng | Spatial Sampling Design for Monitoring Drought, Flood and Earthquake Disaster in China[END_REF].

Referring the bottom left corner of the layer, sampling points are placed according to square grids of different sizes. Using the data of sampling points, a Kriging interpolation is applied to estimate the value of non-sampled cultivated land. Compare the grid sizes and the corresponding MSEs; then a suitable grid size can be chosen for the basic network.

Global monitoring network.

Becasue the center of a grid may be not located inside a map spot of cultivated land, we create a searching circle around this center. Choose the centroid of a map spot, which is closest to the center and inside the searching circle be a candidate point. As these candidate points may raise the estimation variance, a further observation is needed to decide whether to add them or not. The observation should be: apply Kriging interpolation and calculate the MSE after adding each candidate point; add the one that decreases MSE by the least to the basic network and take the points which make MSE decline as candidates for the next calculation. The loop stops when no point reduces the global MSE.

Densification of the network. Adding sampling points in complex areas where natural quality of cultivated land changes intensely can not only increase the accuracy of interpolation but also be quite meaningful because there may be a high chance of quality grades change.

At present, the main approach to densify the network is to use existing points to get weights of non-sampled locations and then add the point with largest weight. A. G. Journel's method was to apply Kriging interpolation at first and get the Kriging variance of every non-sampled point, then set the point with largest Kriging variance as the new one to add [START_REF] Journel | Nonparametric Geostatistics for Risk and Additional Sampling Assessment[END_REF]. Kriging variance is used to evaluate the uncertainty of Kriging interpolation [START_REF] Van Groenigen | Constrained Optimisation of Soil Sampling for Minimisation of the Kriging Variance[END_REF]. But in this study, real values of cultivated land blocks were known so that error of estimation, which shows the deviation of real and estimated value, was more reliable than Kriging variance. Hence, estimation variance was used as the weight in the study. However, points with large estimation variance may stand as outliers; to solve this problem, tessellation polygon was introduced. Tessellation polygon was proposed by Netherland meteorologist A. H. Thiessen to predict average rainfall based on data of weather stations in discrete distribution. He linked all adjacent weather stations as triangles and drew perpendicular bisector of each line to produce tessellation polygon. At last, the rainfall intensity of polygon region was presented by the only weather station inside the polygon. The tessellation polygon has the following features: a) There is only one discrete point inside each polygon. b) All objects inside a polygon are closest to the corresponding discrete point. c) The objects on the line of a polygon have the same distance to the two points on both sides of the line.

When using monitoring points to build tessellation polygon, every polygon is the control area of one point. If MSE of some polygon is relatively high, the corresponding estimation accuracy is low. To some degree, this approach can exclude individual outliers to be added. The details are listed as follows:

Selecting the point of maximum variance in a tessellation polygon with a high MSE. Add it to the network only if it makes the global MSE smaller. The procedure of densification stops until the count of sampling points reaches the specific value or the MSE of tessellation polygon is lower than the threshold.

Till then, the work to set up the monitoring network is complete.

Estimate the population mean

After monitoring network is finished, the population mean of natural quality indices of cultivated land will be calculated by following formula:

       n 1 i x n 1 x i i f f (5)
Where n means the total of cultivated land map spots, which is 5182;   i f x is the population mean;   i f x ˆis the estimation value.

Results and Discussion

Set up network based on Kriging estimation variance

According to the procedure mentioned above, different sizes of square grid were used and searching circle's radius was set as 100m. In principle, the smaller the size of grid is, the larger the count of points is. But in light of discontinuous distribution of cultivated land, the count was affected by the position relation of grids and map spots. So, it is legal that sometimes smaller grids contain fewer points. Table 1 shows the number of points and MSE of each grid size.

Considering both MSE and number of points, we found it reasonable to choose the 4000m grid for three reasons: (a) The MSE of 4000m grid with 39 points was 13130, smaller than that of 3900m, 3800m, 3600m, 3500m, 3300m and 3100m. This meant that the 4000m grid contained fewer points but had higher accuracy. (b) Although the MSEs of 3700m, 3400m, 3200m or sizes below 3000m were smaller, the number of points increased significantly along the way, followed by the augmentation of sampling cost. (c) To leave basic monitoring network for further densification, we made the number of points as small as possible. Even though the searching radius was set up as 100m, from the figure 4(a), we can see that sampling points did not cover every grid. To fix this, we set the radius as 500m and got 11 candidate points, as showed in figure 4

Results Comparison

Compared with classical sampling models.

(a) Simple Random Sampling Table 2 shows twenty groups of results using simple random sampling with 48 points to set up the network. According to the table, the minimal MSE was 13171.69, which was bigger than that of our method. When applying stratified sampling, we divided study area into 4 stratums (from 11 to 14 natural quality grade) and placed 48 sampling points as well. The formula followed [START_REF] Jinfeng | Spatial Sampling and Inference[END_REF] was used to assign the number of samples for each stratum: The results indicated classical sampling models lead to low estimation accuracy of interpolation because of the correlations among samples.

, n C W C W n h h h h h    h h S S N N W h h  (
Compared with traditional square grids. From table 1, the count of points of 3800m and 3900m equaled that of our method, but the estimation accuracy was significantly lower than ours. Besides, until sampling points reached one hundred, traditional grids can provide sufficient accuracy. In summary, the densification of our method is valid.

Compared with adding points without tessellation polygon. If we did not use tessellation polygon as a limit and directly added points of maximum variance, the MSE was 11680.68.Due to the fact that point of maximum variance could be quite different from surrounding values, directly adding may decrease the estimation accuracy and MSE did not monotonically decrease during the procedure of densification. This experience illustrated that accuracy could be improved when we limited the new points to the area of tessellation polygon with relatively high MSE.

Conclusion

By comparing our approach to traditional sampling methods in cost (count of sampling points), surface fitting (measured by MSE) and inference of the population mean, we can conclude that:

(a) Because spatial samples are correlated, traditional sampling methods will reveal their shortages like low efficiency and accuracy. We came up with a method based on Kriging estimation variance, which takes the spatial variation and distribution characteristics into account, and involves interpolation to calculate the value of non-sampled cultivated land.

(b) When densifying the network, points of maximum variance cannot be added directly in case that some outliers appear to decrease accuracy.

Introducing tessellation polygon can deal with the problem to some degree.

(c) The paper reports the way to build the monitoring network for cultivated land quality in Beijing Daxing district. By interpolating the whole area, the estimation accuracy is better than traditional methods and the statistical inference of population mean stays in an accurate level as well.

  consideration are as follows: a) Natural quality index could reflect the natural conditions and agricultural facilities of the region. b) After calculating the natural quality index, the corresponding utilization and economy quality indices will be obtained by introducing the coefficients which are assigned by Chinese provinces [8]. c) Utilization and economy quality indices are associated with input cost and grain output, so they are largely affected by human factors and spatial interoperation cannot be applied directly.

Fig. 1 .

 1 Fig.1. Study area. Daxing district, located in the south of Beijing, has four natural quality grades (from 11 to 14), showing the ladder-like distribution from northeast to southwest. Higher quality land of 14th grade mainly in northeast makes up 7.72%, that of 13th grade 36.19%, that of 12th grade 35.35% and the rest mainly located in middle and south area belongs to 11th grade.

  Data examination is necessary before Kriging interpolation because Kriging is based on the assumption of stationarity, which demands all data values come from distributions that have the same variability. Secondly, certain Kriging methods require the data to be approximately normally distributed. Normal QQ Plot was used to examine the data in the study. In the diagram below, Data points fall close to the 45-degree reference line, approximating normal distribution. Therefore, the application of Kriging interpolation is allowed.

Fig. 2 .

 2 Fig.2. Normal QQ Plot: the quantile values of the standard normal distribution are plotted on the x-axis, and the corresponding quantile values of the dataset are plotted on the y-axis.

Fig. 3 .

 3 Fig.3. Workflow of setting up the network: N could be set to the final count of samples according to the cost (N=5182by default). A threshold could be set and new points are added to n regions whose MSEs are bigger than the threshold (n equals the number of global sampling points by default). Variable m is initially set to 1.

Fig. 4 .

 4 Fig.4. (a) is the basic network of 39 sampling points by 4000m grids. (b) shows the locations of candidate points. (c) is the global network of 42 points. (d) shows the tessellation polygons.

Fig. 5 .

 5 Fig.5. Distribution of monitoring points in Daxing district

Table 1 .

 1 Number of points and MSE of each grid size Grid size Number MSE Grid size Number MSE Grid size Number MSE

	4500	34	17077.8 3600	53	13258.8 2700	88	12288.0
	4400	34	15527.0 3500	59	13288.6 2600	94	10338.6
	4300	36	13763.5 3400	49	12336.4 2500	108 12510.7
	4200	37	13995.1 3300	49	13635.9 2400	102 11123.3
	4100	31	15932.5 3200	67	12913.5 2300	126 10857.5
	4000	39	13130.1 3100	74	14568.0 2200	141 10602.5
	3900	48	13930.6 3000	69	12693.5 2100	154	9505.5
	3800	48	14675.4 2900	74	12402.4 2000	161 10046.4
	3700	51	12413.1 2800	82	11072.8

Table 2 .

 2 Results of simple random sampling

	Times	MSE	Times	MSE	Times	MSE	Times	MSE
	1	13423.81	6	13759.86	11	13171.69	16	14000.15
	2	14775.30	7	13222.51	12	15108.00	17	14467.08
	3	15999.14	8	14419.24	13	15345.55	18	14091.05
	4	14154.64	9	14595.11	14	16037.21	19	13274.75
	5	13180.20	10	18165.41	15	14811.47	20	15200.25
	(b) Stratified Sampling

6 )

 6 Where h means the hth stratum; N stands for total number of cultivated land; N Based on the assigned numbers in table 3, the study used random method to place them. Results of twenty experiences are contained in table 4, with the minimal MSE being 13990.60 and maximal one being 18546.03, far worse than 11238.32 of our method.

h is the number of land in hth stratum; W h denotes the weight of hth stratum; S h means standard deviation of natural qualities of h stratum; And C h is the cost for setting up a sampling point.

Table 3 .

 3 Standard Deviation and number of samples in each stratum

	Natural quality classification	11	12	13	14
	Standard Deviation	19.37	26.57	49.16	19.49
	Number of sampling points	7	13	25	3

Table 4 .

 4 Results of stratified sampling

	Times	MSE	Times	MSE	Times	MSE	Times	MSE
	1	15107.41	6	17697.43	11	14791.12	16	18546.03
	2	13990.60	7	18089.52	12	14441.22	17	17209.28
	3	16501.53	8	15171.63	13	14660.60	18	14002.29
	4	17101.16	9	16482.26	14	14697.30	19	16447.92
	5	14936.49	10	14722.58	15	16422.40	20	16319.17

the information management system for Farmland grading monitoring (201011006-5).


Zhimei Li

Ningbo Wang

Jiajiao Wu

Jay Richard Stauffer

Zhihong Li

The Potential Geographical Distribution of Bactrocera cucurbitae (Diptera: Tephritidae) in China Based on Eclosion Rate Model and ArcGIS

Keywords: Bactrocera cucurbitae, potential geographical distribution, ArcGIS, eclosion rate model, plant quarantine

The melon fruit fly, Bactrocera cucurbitae Coquillett (Diptera: Tephritidae), is one of the important insect pests of fruits and vegetables. In order to monitor and control it effectively, it is necessary to know the potential geographical distribution of this pest. The ER (Eclosion rate) model was constructed from empirical biological data, and analyzed with ArcGIS. Based on the soil temperature and moisture data of Chinese meteorological stations, the potential geographical distribution of B. cucurbitae from January to December in China was predicted. Six categories were used to describe different levels of suitability for B. cucurbitae in China. The potential geographical distribution and suitable levels for every month in China were obtained and showed that almost all locations were suitable from May to September. Further analysis showed that monitoring measures should be taken in Guangdong, Guangxi, Yunnan, and Hainan provinces throughout the year.

INTRODUCTION

The melon fruit fly, Bactrocera cucurbitae Coquillett (Diptera: Tephritidae), originated from India [START_REF] Dhillon | The melon fruit fly, Bactrocera cucurbitae: A review of its biology and management[END_REF], spread throughout Southeast Asia in the 1990s, and then spread to other regions. At present, B. cucurbitae is found in some areas of Asia, Africa, and the Pacific Islands, including more than 30 tropical and subtropical countries and regions [START_REF]Animal and Plant Quarantine Administration of the People's Republic of China and Plant Quarantine Institute of Ministry of Agriculture: A Selection of the entry plant quarantine pests of China[END_REF][START_REF] Eppo | Distribution Map of Quarantine Pests for Europe Bactrocera cucurbitae[END_REF]. In China, B. cucurbitae was mainly distributed in Guangdong, Guangxi, Hainan, Fujian, Yunnan (the geographic distribution in Yunnan can be divided into year-round distribution and seasonal distribution. The year-round distribution area was located at the South of latitude 24°), Taiwan, and Hong Kong. In addition, B. cucurbitae was found in Sichuan, Chongqing, Guizhou, Zhejiang, Jiangsu, Hunan, and other places in recent years [START_REF] Wang | Overview about important fruit fly pests of fruits and vegetables[END_REF][START_REF] Wang | Diptera Tephritidae insects of East Asia[END_REF][START_REF] Zhou | Biology and other comprehensive prevention and control of B. dorsalis for mangos[END_REF][START_REF] Jiang | Study on the biology of Bactrocera cucurbitae[END_REF][START_REF] Jiang | Pest Risk Analysis of Bactrocera cucurbitae[END_REF][START_REF] Jiang | Life tables of the melon fly laboratory population reared on various host fruits[END_REF][START_REF] Zhu | Genetic relationships among four Bactrocera cucurbitae geographic populations in Yunnan Province[END_REF]. Bactrocera cucurbitae has been reported to attack more than 100 species of fruits and vegetables, and has been listed as an important quarantine pest by many countries including China [START_REF] Chen | On the seasonal population dynamics of melon and oriental fruit flies and pumpkin fly in Guangzhou area[END_REF][START_REF] Li | Damage and Control of Bactrocera cucurbitae[END_REF].

Bactrocera cucurbitae often caused losses of 30% to 90%. According to the reports, B. cucurbitae had caused a reduction in production of watermelon of India, bitter gourd of Nauru, and towel gourd and pumpkin of Solomon by 28.55%, 95%, 90% and 87%, respectively [START_REF] Dhillon | The melon fruit fly, Bactrocera cucurbitae: A review of its biology and management[END_REF]. Japan had spent 20.4 billion yen to eradicate B. cucurbitae in Okinawa islands [START_REF] Koyama | Eradication of the melon fly, Bactrocera cucurbitae, in Japan: importance of behavior, ecology, genetics, and evolution[END_REF].

Population fluctuation of B. cucurbitae not only has been correlated with the temperature, but also with the relative humidity and rainfall [START_REF] Zaman | Assessment of the male population of the fruit flies through kairomone baited traps and the association of the abundance levels with the environmental factors[END_REF].Temperature had a significant impact on B. cucurbitae [START_REF] Manzar | Population fluctuation of fruit flies, Bactrocera spp. infesting bitter gourd in Central Uttar Pradesh[END_REF], under laboratory conditions; when the temperature was below 14℃, all instars of B. cucurbitae did not develop normally [START_REF] Keck | Effect of temperature on development and activity of the melon fly[END_REF][START_REF] Liang | Fruit fly and its Control[END_REF]. The low temperature of developmental threshold for pupae was 11.14℃, developmental temperature range was 11-36℃ [START_REF] Liang | Fruit fly and its Control[END_REF], and the optimum developmental temperature range was 22-30℃ [START_REF] Yuan | Influence of Temperature on Development, Survival and Fecundity of Bactrocera dorsalis Hendel[END_REF]. When the temperature was greater than 32.2℃, the mortality rate to increased, and under the condition of 34℃ pupae mortality rate was 95.27%±0.35% [START_REF] Yuan | Influence of Temperature on Development, Survival and Fecundity of Bactrocera dorsalis Hendel[END_REF]. Soil moisture also affected the pupae greatly. In India, the report showed two peak period of melon flies attacking the fruits --the one is in June with the relative humidity 72.0% -84.0 % and rainfall 53.3-61.5 mm, the other is the second and third weeks in July with the relative humidity 85.0% and rainfall 114.1-247.5 mm [START_REF] Ingoley | Seasonal occurrence of fruit fly, Bactrocera cucurbitae Coq., on cucumber in midhills of Himachal Pradesh[END_REF]. Mature larvae of B. cucurbitae pupated in the soil of 0.5-15 cm depth [START_REF] Jackson | Depth of pupation in four species of fruit flies (Diptera: Tephritidae) in sand with and without moisture[END_REF][START_REF] Pandey | Studies on movement of Dacus cucurbitae maggot for pupation[END_REF], and it is favorable for the mature larvae to pupate with higher eclosion rate when the soil moisture is less than 25%. When the soil moisture is greater than 30%, both the pupation rate and eclosion rate are lower [START_REF] Yuan | Influence of Temperature on Development, Survival and Fecundity of Bactrocera dorsalis Hendel[END_REF].

The potential geographical distribution study of fruit flies started with, climate maps, ecological niche modeling software, technical methods, etc. and involved several species (e.g., Ceratitis capitata Wiedemann, B. dorsalis (Hendel), B. (Daculus) oleae). In 1924, [START_REF] Cook | The distribution of the alfalfa weevil Phytonomus Posticus Gyll. A study in physical ecology[END_REF] used climate maps predict the pest potential geographical distribution, and then this technology was further developed [START_REF] Cook | The distribution of the alfalfa weevil Phytonomus Posticus Gyll. A study in physical ecology[END_REF]. [START_REF] Messenger | Bioclimatology and prediction of population trends[END_REF] The highly suitable areas of B. cucurbitae Coquillett in China included Guangdong, Guangxi, Hainan, southern of Fujian, southern of Yunnan, western of Taiwan, and the Sichuan Basin, while moderate to low suitable areas include parts of Jiangxi, Hunan, Guizhou, Chongqing, Shanghai and Sichuan, Yunnan, Fujian, Zhejiang, Jiangsu, Anhui, Hubei, Shaanxi, Henan, and Gansu [START_REF] Kong | Effects of climatic factors on the growth and geographical distribution of Oriental fruit fly, Bactrocera dorsalis[END_REF].

The purpose of this study was to establish the ER (eclosion rate) model of B.

cucurbitae, and analyze the potential geographical distribution and suitability levels using the ER model and displaying distribution patterns using soil temperature and soil moisture of past years in China by ArcGIS.

MATERIALS AND METHODS

The B. cucurbitae samples collected from Huangpu in Guangdong province were selected for the eclosion experiment. Eggs of B. cucurbitae were obtained from adults that had been reared for four generations on an artificial diet. Mature larvae (6 days after egg hatch) were placed in moist sand (75% relative humidity) at 29℃ for pupation. All 6480 pupas were gathered after 24h under and held at 25℃. The ER model, based on soil temperature and relative humidity was obtained from a crossover design experiment conducted at the plant quarantine laboratory of Guangdong Entry-Exit Inspection and Quarantine Bureau. The design specified six soil temperature grades: 9℃, 14℃, 19℃, 24℃, 29℃, and 34℃. Six relative humidity grades included 0%, 20%, 40%, 60%, 80%, and 100%. Experiments consisted of 36 treatments and 3 replications during about 25 days and every box had 60 pupas. All pupas were placed 2 cm under the soil and held in artificial climate boxes. Water was added as needed [START_REF] Wu | Insect ecology experiments[END_REF].The ER of B. cucurbitae was obtained and the ER model was derived by stepwise regression (SPSS 13.0) Geographic distribution and suitability was calculated for each month of the year for all of China and displayed as maps by ArcGIS. Soil temperature and moisture data were limited, but at data for a minimum of 10 days per month were obtained. Since values were obtained from different locations in different months, the full data set is presented in Table 4. The suitability of B. cucurbitae was analyzed on the basis of these data and ER model results. The ER of B. cucurbitae for every location in China was obtained. For each location, suitability was categorized into 6 levels; negligible (where the B. cucurbitae is unable to occur and survive) (ER=0), extremely low (0< ER≤0.1), low(0.1<ER≤0.2), moderate (0.2<ER≤0.3), high (0.3<ER≤0.4) and extremely high (0.4<ER≤1). Suitability maps were plotted for each month using the inverse distance weight (IDW) raster interpolation.

RESULTS

ER model

SPSS was used to build the modelas it is convenient, fast, and reliable. The binary model was not significant (p>0.05) and the binary cubic model is too complicated to explain, so we chose binary quadratic regression model as the eclosion rate model of B. cucurbitae. After XY was dropped, the model chosen for eclosion rate was as follows:

22 Z=-0.0041X -0.00005Y +0.17976X+0.01044Y-1.618 Z is the ER (eclosion rate) of B. cucurbitae X is the soil temperature Y is the soil moisture Table 1 shows the regression model summary, and four constants are available including XX, X, Y and YY. Table 2 shows the analysis of variance of regression.

Stepwise regression ordered the independent variables according to their explanatory power. Analysis showed that XY did not contribute to the explanatory power of the model significantly (p>0.05), so it was dropped. By regression analysis, after XY was dropped, t=0.107, P>0.01(Table 3). The model was reliable because of R 2 =0.706, F=48.991>F0.05 and P<0.01. The regression analysis of the ER showed that P was significant on 0.01 level.

Potential geographical distribution of B. cucurbitae in China

After the suitable range and suitable level analysis of B. cucurbitae in China monthly the results showed that the percent of suitable locations was about 11% in January and February, then it increased to 37.14% quickly in March, and it was 76.19% in April.

Almost all locations were suitable from May to September, but it was 61.76% in October, and then decreased to 19% in November, and it was 12.28% in December.

For the most part, it was basically the same in December, January and February, and there was a turning point in March and April with the suitable locations increasing.

May to September was the most suitable period, and there was another turning point in October and November with the suitable locations decreasing. Finally, it accomplished an alternating cycle of the four seasons. Combined with the change of the different suitable level locations, May and September took the largest proportion of extremely high and high suitable level locations, followed by them was April, June, August and October, and then it was July, November, December, February and January. The percentage of extremely high level locations was the lowest and the percent of extremely low locations was the highest in January. Ba county in March when monthly average soil temperature was 10.89°C, the monthly average soil relative humidity was 16.90% (Fig. 1).

As a turning point, we can clearly see the division of the potential geographical

Conclusions

The results indicated that soil temperature and soil moisture affected the potential geographical distribution of B. cucurbitae in China and suitable areas changed with

DISCUSSION

Monitoring and trapping periods

This study indicated four monitoring and trapping periods for fruit fly eclosion in China. If continuous two months are not suitable, monitoring should be suspended.

Monitoring measures should be taken throughout the year in Guangdong, Guangxi, Yunnan, and Hainan. The second period was from March 1 to November 30 and included Hubei, Hunan, Sichuan, Chongqing, Guizhou, Jiangxi, and Fujian. The third period, from March 1 to October 31 included Shandong, Jiangsu, Shanghai, Zhejiang, Anhui, Tibet. and Taiwan. The last period was from April 1 to October 31, and monitoring areas included Liaoning, Heilongjiang, Jilin, Inner Mongolia, Qinghai, Beijing, Tianjin, Hebei, Henan, Shaanxi, Shanxi, Gansu, Ningxia, and Xinjiang.

In summary, the two important natural factors (soil temperature and moisture) determined the potential geographical distribution of B. cucurbitae over 12 months in China were reported, and the ER model can be used to predict where and when these conditions are suitable for eclosion. This model also can be applied to other countries and will be improved with the development of the future research.

Potential geographical distribution of B. cucurbitae in China

  used artificial climate chambers to simulate a dozen typical weather conditions of the United States, to study the growth and development of the C. capitata Wiedemann, B. dorsalis (Hendel) and B. cucurbitae under different climatic conditions. Combined with climate analysis, the suitable geographical distribution of the three fruit fly species in the United States was described [23]. Yonow applied CLIMEX to study the potential geographical distribution of the B. (B.) tryoni in Australia [24]. The potential geographical distribution study of fruit flies started later in China, but there are many papers in recent years, including the Multiple Indicator fuzzy comprehensive evaluation technology, agro-climatic similarity distance technology, geographic information systems ( GIS ) technology and ecological niche modeling software, among others. Many fruit flies were involved such as B. dorsalis (Hendel), B. cucurbitae Coquillett, Anastrepha ludens Loew and Rhagoletis pomonella Walsh. Zhou (2005) built the lethal temperature model and the effective accumulated temperature model of B. cucurbitae Coquillett in 2005, and used 30 years' climate data of the 670 meteorological stations in China to run the model. The prediction results showed that the melon fly could occur in 36.27% of China, the northern boundary of the distribution is about 31°±2°N. It could produce 2-10 generations/year, and mainly produced 4-7 generations/year [25]. Kong et al. (2008) analyzed the potential geographical distribution of B. dorsalis (Hendel) and B. cucurbitae Coquillett in China and world by combining the CLIMEX and DIVA-GIS in 2008.

  distribution of B. cucurbitae in October. Extremely high and high suitable area concentrated in the southern provinces of China, in which we must pay attention to the Nuomuhong of Qinghai and Tianjin. Yili of Xinjiang (latitude 43.57°, longitude 81.20°) was the northern boundary of the suitable regions in October. Monthly average soil temperature was 10.93ºC , and the monthly average soil relative humidity was 17.01% (Fig. 2).

Fig. 1 .Fig. 2 .

 12 Fig. 1. The potential geographical distribution of B. cucurbitae in March

Table 1 .

 1 Regression model summaries

				Adjusted R	Std. Error of the
	Model	R	R Square
				Square	Estimate
	1	0.840(a)	0.706	0.692	0.206
	a Predictors:(Constant), xy, xx, yy, y, x

Table 2 .

 2 Analysis of variance of regression

	Model		Sum of squares	Df	Mean square	F	Sig.
	1	Regression	10.396	5	2.079	48.991 1.31E-25
		Residual	4.329	102	0.042

Table 3 .

 3 Excluded variables of regression model

	Model		Unstandardized Coefficients	Standardized Coefficients	t	Sig.
			B	Std. Error	Beta
	(Cons
	tant)	-1.618	0.158		-10.218 2.69E-17
	x		6.112	0.485	4.157	12.593	1.68E-22
	y		1.004	0.253	0.928	3.961	0.0001
	xx		-4.743	0.367	-4.138	-12.906 3.59E-23
	yy		-0.539	0.199	-0.519	-2.713	0.0078
	xy		0.025	0.231	0.018	0.107	0.9152
	a	Dependent Variable: z

  The percentage of extremely high and high level locations increased linearly from January to May, and then it decreased from May to July. It increased from July to September, and then decreased from September to December. It was the lowest in July because the soil temperature exceeded the eclosion threshold temperature of B. cucurbitae. The distribution of B. cucurbitae diffused towards north in March. Guangdong, Guangxi, Yunnan and Sichuan came to be the extremely high or high suitable level locations. We must pay attention that Shannan prefecture of Tibet, Juxian and Laiyang of Shandong came to be the extremely high or high suitable level locations in March. Ba county of Hebei Province (latitude 39.07°, longitude 116.23°) was the northern boundary of the suitable regions in March. The eclosion of B. cucurbitae occurred at

ACKNOWLEDGMENTS

The authors would like to thank Liang Guangqin, Mo Renhao, Hu Xuenan, Liu Haijun, Li Chunyuan and all other officers of Guangdong Inspection and Quarantine Technology Center to support the experiment. Thanks for all members of the Plant Quarantine and Invasion Biology Laboratory of China Agricultural University (CAUPQL). This study was supported by the Program of Ministry of Agriculture, China (No. 2012-Z15).

  * This study was supported by the Program of Ministry of Agriculture, China (No. 2012-Z15).

Almost all suitability analysis methods of fruit flies include the main natural factors that affect the occurrence of fruit flies, and temperature and humidity often is the key to determine fruit fly outbreaks. Under the situation that meteorological data are relatively limited, how to make full use of these data is particularly important; however, we often encounter data substitution or conversion problems (e.g., temperature and soil temperature data conversion), as well as the relationships among atmospheric humidity, rainfall, and soil moisture, which are very important to the adaptability analysis research of fruit flies.

In addition, the mode of 'biological experiments + biological model +GIS analysis' was suitable for suitability analysis of invasive fruit flies. To date, pest suitability analysis techniques use models combined with software analysis as its main analysis methods; thus, models and software for medical and military fields can be used to analyze the pest suitability. Once related environmental meteorological data are available, we can analyze the suitability of any instars of the invasive species year by year, month by month, even ten days by ten days and day by day, and a detailed study can increase the accuracy of suitability analysis and forecast. Detailed and completed information obtained from the literatures and biological experiments is very important for the potential geographical distribution research. We also need to explore how to combine the results of quantitative analysis and graphics, as the analysis of the species turns from qualitative methods towards quantitative methods.

Management countermeasures

Quarantine measures should be strengthened domestically as well as at the borders (inspection and quarantine). The critical tasks of border inspection and quarantine are port inspection and quarantine from May to September, because many areas of China are suitable for B. cucurbitae during this time. Especially from imports of these continents in which the B. cucurbitae distributes: Asia (Oman, Bangladesh, Brunei, Cambodia, India, Indonesia, Japan's Ryukyu Islands, Laos, Malaysia, Myanmar, Nepal, Pakistan, Philippines, Sri Lanka, Thailand, Vietnam), Africa (Egypt, Kenya, Tanzania, Mauritius, some areas of Reunion), and Solomon Islands, Mariana Islands, Papua Island, Hawaii of USA and other parts of the Pacific Islands [START_REF] Dhillon | The melon fruit fly, Bactrocera cucurbitae: A review of its biology and management[END_REF][START_REF] Eppo | Distribution Map of Quarantine Pests for Europe Bactrocera cucurbitae[END_REF][START_REF] Stonehouse | The Melon Fly Bactrocera cucurbitae in the Republic of Seychelles[END_REF]. In addition, domestic inspection and quarantine should monitor the transport of fruits and vegetables from south to north in China, such as Guangdong, Guangxi, Hainan, Fujian, Yunnan, Taiwan, Hong Kong, Sichuan, Chongqing, Guizhou, Zhejiang, Jiangsu, Hunan, and other places [START_REF] Wang | Overview about important fruit fly pests of fruits and vegetables[END_REF][START_REF] Wang | Diptera Tephritidae insects of East Asia[END_REF][START_REF] Zhou | Biology and other comprehensive prevention and control of B. dorsalis for mangos[END_REF][START_REF] Jiang | Study on the biology of Bactrocera cucurbitae[END_REF][START_REF] Jiang | Pest Risk Analysis of Bactrocera cucurbitae[END_REF][START_REF] Jiang | Life tables of the melon fly laboratory population reared on various host fruits[END_REF][START_REF] Zhu | Genetic relationships among four Bactrocera cucurbitae geographic populations in Yunnan Province[END_REF]. Quarantine measures ought to be reinforced to prevent the infestations of fruits and vegetables.


Xiongbing Tu

Zhihong Li

Zehua Zhang

Zhigang Wu

Wenlong Ni

Liao Fu

Yasen • Shali

The potential geographical distribution of Locusta migratoria tibetensis Chen (Orthoptera: Acrididae) in Qinghai-Tibet Plateau

Keywords: Locusta migratoria tibetensis Chen, CLIMEX, potential geographic distribution, ArcGIS, climate change

Locusta migratoria tibetensis Chen (Orthoptera: Acrididae) is one of the most harmful species of Acrididae. It causes extensive damage to gramineous crops. The climate mapping program, CLIMEX 3.0, and the GIS software, ArcGIS 9.3, were used to model the current and future potential geographical distribution of L. m. tibetensis. The model predicts that, under current climatic conditions, L. m. tibetensis will be able to establish itself throughout much of Qinghai-Tibet Plateau in China, including Qinghai, Tibet and Sichuan province. A total of 94.23% of Qinghai-Tibet Plateau's land mass or 2.22×10 6 km 2 , is climatically suitable under current climatic conditions. In this study, the potential distribution of L. m. tibetensis under future climate scenarios was projected and quantified climatically suitable areas. In Qinghai-Tibet Plateau, the total suitable climate land mass was increased by 0.16×10 6 km 2 or 0.70% until 2100. However, compared to 2020 and 2050, the suitable land mass was decreased by 0.47×10 6 km 2 and 0.17×10 6 km 2 . The most optimal area would increase 0.16 ×10 6 km 2 until 2100, the most of Qinghai-Tibet Plateau would be expected to offer optimal climatic conditions

Introduction

Locusta migratoria tibetensis Chen is one of the most harmful pests and causes large amounts of damage to gramineous crops. Since 1990s, the occurrence of L. m. tibetensis has showed a rising trend. Especially in 2003 to 2006, there were 1000 to 3000 individuals per m 2 along the valley areas of Jinsha River, Yarlung Zangbo River and Yalong River, which had threaten the production of highland barley and grass [START_REF] Chen | Historical evidence for population dynamics of Tibetan migratory locust and the forecast of its outbreak[END_REF][START_REF] Feng | The change of potential distribution of Locusta migratoria tibetensis Chen (Orthoptera: Acrididae) from 1961 to 2005[END_REF] .

The distribution area of L. m. tibetensis is only in China at present, including northern of Qinghai, Tibet and west of Sichuan province [START_REF] Chen | A new subspecies of Locusta migratoria-Tibetan migratory locust (Locusta migratoria tibetensis subsp. n.)[END_REF] . Studies showed that there was a significant correlation between increasing of potential distribution area and global warming. The annual surface temperature increased at the speed of 0.0301℃, while the potential distribution area increased 504.38 km 2 per year. This result indicated that a little increment in temperature would bring a great extension in distribution of L. m. tibetensis [START_REF] Feng | The change of potential distribution of Locusta migratoria tibetensis Chen (Orthoptera: Acrididae) from 1961 to 2005[END_REF] .

The CLIMEX model has been used to describe the potential distribution of some pests, such as Ceratitis capitata (Wiedemann) [START_REF] Worner | Ecoclimatic assessment of potential establishment of exotic pests[END_REF][START_REF] Vera | Potential geographical distribution of the Mediterranean fruit fly, Ceratitis capitata (Diptera:Tephritidae), with emphasis on Argentina and Australia[END_REF] , Bactrocera tryoni(Froggatt) [START_REF] Yonow | The geographical distribution of the Queensland fruit fly, Bactrocera tryoni, in relation to climate[END_REF][START_REF] Sutherst | The vulnerability of Australian horticulture to the Queensland fruit fly, Bactrocera tryoni, under climate change[END_REF] ,

Carpomya vesuviana Costa [START_REF] Lv | Potential geographic distribution of Ber fruit fly, Carpomya vesuviana Costa, in China[END_REF] , B. tsuneonis(Miyake) [START_REF] Wang | The potential geographic distribution of Bactrocera tsuneonis (Diptera: Tephritidae)[END_REF] , B. scutellata (Hendel) [START_REF] Ni | Potential geographic distribution of Bactrocera scutellata (Hendel), in China[END_REF] , Toxotrypana curvicauda Gerstaecker [START_REF] Qu | Potential geographic distribution of papaya fruit fly, in China[END_REF] , B. zonata(Saunders) [START_REF] Ni | Including climate change in pest risk assessment: the peach fruit fly, Bactrocera zonata (Diptera: Tephritidae)[END_REF] , Melanoplus sanguinipes (Fabricius) [START_REF] Olfert | Application of general circulation models to assess the potential impact of climate change on potential distribution and relative abundance of Melanoplus sanguinipes (Fabricius) (Orthoptera: Acrididae) in north American[END_REF] . This model has also previously been used to predict the effects of climate change on species' potential distributions using either a regional global climate model [START_REF] Kriticos | The role of modelling in weed management[END_REF] or synthetic climates [START_REF] Kriticos | Climate change and biotic invasions: a case history of a tropical woody vine[END_REF][START_REF] Kriticos | Climate change and the potential distribution of an invasive alien plant: Acacia nilotica ssp. indica in Australia[END_REF] . Species examined using this approach include B. dorsalis (Handel) [START_REF] Stephens | The current and future potential geographical distribution of the oriental fruit fly, Bactrocera dorsalis (Diptera: Tephritidae)[END_REF] , Melaleuca quinquenervia [START_REF] Watt | The current and future potential distribution of Melaleuca quinquenervia[END_REF] and Melanoplus sanguinipes (Fabricius) [START_REF] Olfert | Application of general circulation models to assess the potential impact of climate change on potential distribution and relative abundance of Melanoplus sanguinipes (Fabricius) (Orthoptera: Acrididae) in north American[END_REF] . To know the potential distribution of the L. m. tibetensis under climate change assumptions due to potential climate change impact, the CLIMEX 3.0 model was used to assess the response of L. m. tibetensis to climate and to make predictions for the 2100s.

Materials and methods

Overview of the CLIMEX model

CLIMEX is a dynamic model [START_REF] Sutherst | CLIMEX Version 3: User's Guide[END_REF] that integrates the weekly responses of a population to climate using a series of annual indices. The growth and stress indices are calculated weekly and then combined into an overall annual index of climatic suitability, the Ecoclimatic Index (EI). This index is expressed on a scale of 0 to 100. Some suggested guidelines are provided [START_REF] Sutherst | CLIMEX v2 CD and User's Guide[END_REF][START_REF] Sutherst | A climate model of the red imported fire ant, Solenopsis invicta Buren (Hymenoptera:Formicidae): implications for invasion of new regions, particularly Oceania[END_REF] based on the EI. An EI of 0.00 indicates that the climate is unsuitable for the species and that the species cannot persist in an area under average environmental conditions. An EI of 0.01-9.99 indicates marginal conditions, an EI of 10.00-19.99 indicates suitable conditions, and an EI of 20+ indicates optimal conditions. An EI of 100 indicates that conditions are perfect throughout the year. However, few environments are sufficiently stable to provide perfect habitat all year [START_REF] Stephens | The current and future potential geographical distribution of the oriental fruit fly, Bactrocera dorsalis (Diptera: Tephritidae)[END_REF] .

To visualize the results, the CLIMEX output was 'loose-coupled' [START_REF] Kriticos | Climate change and the potential distribution of an invasive alien plant: Acacia nilotica ssp. indica in Australia[END_REF] to a geographical information system (ArcGIS). This system was purchased from ESRI to create thematic maps. ArcGIS was also used to make projections of the land areas that would have a suitable climate for the pest species.

Meteorological databases and climate change

Two climate databases were used in this modeling exercise. By default, CLIMEX uses 30-year averages of climate data to estimate climatic suitability. The CLIMEX standard meteorological dataset was first used to create an initial fit. This dataset accompanied CLIMEX version 3.0 and consists of 30-year averages from 1961 to 1990 for an irregularly spaced set of 2500 climate stations. Subsequently, a regular girded dataset of the normal climate for the same period TYN SC 2.0 (http://www.cru.uea.ac.uk/~timm/grid/TYN_SC_2_0.html) [START_REF] Kriticos | Predicting the potential geographic distribution of weeds in[END_REF] was used to fine-tune the parameter fit. The dataset consisted of 67,420 points spaced on a 0.5° latitude×0.5° longitude regular grid for significant land areas worldwide.

The future average annual temperature will increase 1.3-2.1℃, 2.3-3.3℃, 3.9-6℃ until 2020, 2050, 2100 respectively. While the future annual mean precipitation will increase 2-3%, 5-7%, 11-17% until 2020, 2050, 2100 respectively in China [START_REF] Stephens | The current and future potential geographical distribution of the oriental fruit fly, Bactrocera dorsalis (Diptera: Tephritidae)[END_REF][START_REF] Qin | Climate and environment evolution in China[END_REF] . In this study, we use the average data of temperature and annual mean precipitation for the potential geographical distribution analysis in the future.

Fitting CLIMEX parameters

To fit the CLIMEX model for L. m. tibetensis, the parameters were manually and iteratively adjusted until the simulated geographical distribution, as estimated by the EI values, coincided with the species' known native distribution. Parameters used in the CLIMEX model are presented in table 1.

Temperature index

In this paper, the minimum temperature for development (DV0) was set at 14.6°C (Li et al., 2007). The lower and upper optimum temperatures for L. m. tibetensis populations were set at 15-28°C [START_REF] Li | The biological characteristics of Locusta migratoria tibetensis (Chen)[END_REF] . So the lower and upper temperature optima (DV1 and DV2) were set at 15 and 28°C, respectively [START_REF] Li | The biological characteristics of Locusta migratoria tibetensis (Chen)[END_REF][START_REF] Jia | The research about the life habit, classes and biological features of locust at Ali Prefecture of TAR[END_REF] . As a result of the species is established in Tibet, and the day average temperature is no more than 32°C (http://cdc.cma.gov.cn/), so the upper threshold temperature (DV3) was set at 32°C.

Moisture index

In commonly, Locust like desiccated environment for embryonic development stage. They survive well when the soil moisture in May at 0.03-0.3 [START_REF] Guo | The biology of the migratory locusts in China[END_REF] . However, we have found that the soil moisture can be reach to 0.40 in October, this would be related with there was a long snow-covered time in Qinghai-Tibet Plateau. So in this paper, the lower soil moisture limit for development (SM0) was set to 0.02, the lower and upper limits for optimal growth (SM1 and SM2) were set to 0.03 and 0.5 respectively [START_REF] Jia | The research about the life habit, classes and biological features of locust at Ali Prefecture of TAR[END_REF][START_REF] Chen | The main locust and ecological management of locust plague in China[END_REF] . According to one year soil moisture changes of Tibet, the upper soil moisture limit for development, SM3, was set to 0.7 (http://cdc.cma.gov.cn/).

Diapause parameter

Diapause is an important biology characteristic for locust to endure bad environment. In this paper, we take Melanoplus sanguinipes (Fabricius) for setting diapause parameter of L. m. tibetensis [START_REF] Olfert | Application of general circulation models to assess the potential impact of climate change on potential distribution and relative abundance of Melanoplus sanguinipes (Fabricius) (Orthoptera: Acrididae) in north American[END_REF] . Compared to Melanoplus sanguinipes, the diapause day length and development days of L. m. tibetensis are very close to Melanoplus sanguinipes (http://cdc.cma.gov.cn/), also L. m. tibetensis is winter diapauses [START_REF] Guo | The biology of the migratory locusts in China[END_REF] . However, there was no related report about diapause induction and termination temperatures.

Stress parameter

All stress parameters were adjusted by the real distributions in China, to ensure that the predictions of areas suitable for the insect were consistent with the insect's current absence.

CLIMEX model

The land areas projected to have a climate suitable under current and future conditions for L. m. tibetensis (EI>0) were quantified for Qinghai-Tibet Plateau (table 2).

Table 2. Projected area that is climatically suitable (EI>0) for L. m. tibetensis under the current climate, expressed as an area (10 6 km 2 ) and as a percentage of the total land of Qinghai-Tibet Plateau, and the percentage changes in Qinghai-Tibet Plateau for 2020, 2050 and 2100 in comparison with those obtained under the current climate.

The climate model was run for the future climate changes in China [START_REF] Qin | Climate and environment evolution in China[END_REF] .

Area

Results

The potential distribution of L. m. tibetensis under current climatic conditions

The projections using the model suggest that under current climatic conditions, L.

m. tibetensis can potentially establish itself to the Qinghai-Tibet Plateau (fig. 1,a).

The results of the analysis indicate that Qinghai-Tibet Plateau would be expected to offer optimal climatic conditions for the possible spread and establishment of the pest.

After CLIMEX modeling, the extent of the land area that is climatically suitable for L.

m. tibetensis under the current climate is quantified by ArcGIS in table 2. A total of 94.23% of Qinghai-Tibet Plateau's land mass or 2.22×10 6 km 2 , is climatically suitable.

The potential distribution of L. m. tibetensis under future climate scenarios

In this study, the potential distribution of L. m. tibetensis under future climate scenarios was projected and quantified climatically suitable areas. In Qinghai-Tibet Plateau, the total suitable climate land mass was increased by 0.16×10 6 km 2 or 0.70% until 2100 (table 2). However, compared to 2020 (fig. 1,b) and 2050 (fig. 1,c), the suitable land mass was decreased by 0.47×10 6 km 2 and 0.17×10 6 km 2 (table 2). The most optimal distribution area (EI: 20.00+) would increase 0.16 ×10 6 km 2 until 2100 (fig. 1,d).

Discussion

L. m. tibetensis as the independent subspecies which has been nominated by Chen [START_REF] Chen | A new subspecies of Locusta migratoria-Tibetan migratory locust (Locusta migratoria tibetensis subsp. n.)[END_REF] . The taxonomy of L. m. tibetensis has been studied recently; Zhang et al. [START_REF] Zhang | Some unorthodox views on the classification and evolution of the migratory locusts in China prompted by molecular population genetic study[END_REF][START_REF] Zhang | Unexpected relationships of substructured populations in Chinese Locusta migratoria[END_REF] had illustrated the potential importance of phylogeographical history in shaping the divergence and distribution patterns of L. m. tibetensis and other subspecies, he thought that L. m. tibetensis had independent genetic evidence for strong population substructure, as a result of geographical blocking for there was no gene flow between the inner and outside of Qinghai-Tibet Plateau. This is the reason for the current distribution of L. m. tibetensis, and we only predicted the potential distribution in Qinghai-Tibet Plateau in this paper.

The results of this study indicate that L. m. tibetensis is expected to be able to establish itself under temperate zones. Such conditions are to be found in most of Qinghai-Tibet Plateau. Moreover, the results suggested that as climate changes, the potential distribution of L. m. tibetensis would expand into cold regions, as it would expand into northward of the Plateau.

We have known that L. m. tibetensis is only in Qinghai-Tibet Plateau at present [START_REF] Chen | A new subspecies of Locusta migratoria-Tibetan migratory locust (Locusta migratoria tibetensis subsp. n.)[END_REF] . Studies showed that a little increment in temperature would bring great extension in distribution of L. m. tibetensis [START_REF] Feng | The change of potential distribution of Locusta migratoria tibetensis Chen (Orthoptera: Acrididae) from 1961 to 2005[END_REF] . This study presented an example of influence of global warming on plateau ecology. So how climatic changes influence the potential distribution of insect in Qinghai-Tibet Plateau will need more attentions [START_REF] Huang | Studies on the acridoid fauna of Xizang and it's evolution[END_REF][START_REF] Chen | A new record of the Desert locust, Schistocerca Gregaria (Forskal) in Tibet, China[END_REF] .

Potential distributions cannot be predicted based on climate alone. There is also a need to consider dispersal and species interactions, such as host availability. For example, mangoes, Mangifera indica Linn, are distributed in over 70 countries, more than 90% are grown in Asia (India, Pakistan, Bangladesh, Myanmar, Malaysia, China, Thailand, etc.), while Tanzania, Zaire, Brazil, Mexico, the United States, etc., also cultivate mangoes [START_REF] Ye | World mango production and trade[END_REF] . Competition and the effects of natural enemies [START_REF] Baker | The role of climate mapping in predicting the potential geographical distribution of non-indigenous pests under current and future climates[END_REF] should also be considerations. Moreover, the limiting factors considered to affect the geographical distributions of pests should include soil type, geographical features, natural and geographical barriers and human activities. For example, the host plant B. tsuneonis limited citrus [START_REF] Wang | The potential geographic distribution of Bactrocera tsuneonis (Diptera: Tephritidae)[END_REF] , so the prediction of the potential distribution of B. tsuneonis in China considered the factor of host distribution. Additionally, the irrigation and land use map were also taken into consideration for Carpomya vesuviana [START_REF] Lv | Potential geographic distribution of Ber fruit fly, Carpomya vesuviana Costa, in China[END_REF] and Rhagoletis pomonella [START_REF] Geng | Potential geographical distribution of Rhagoletis pomonella (Diptera: Tephritidae) in China[END_REF] . The land use map removed the desert, river and lake, and urban land use type, etc. Therefore, the predictions based on the model and the results discussed above surely have some limitations; these additional factors should also be considered, and more detailed risk assessments should be useful for obtaining more scientific predictions.

Table 1 . 4

 14 Parameters used in the CLIMEX model for the L. m. tibetensis Parameter Value DV0 Lower threshold temperature DV1 Lower optimum temperature DV2 Upper optimum temperature DV3 Upper threshold temperature SM0 Lower threshold of soil moisture SM1 Lower limit of optimum soil moisture SM2 Upper limit of optimum soil moisture SM3 Upper threshold of soil moisture DPD0 Diapause induction day length DPT0 Diapause induction temperature (average weekly minimum) DPD1 Diapause termination temperature (average weekly minimum) DPD Diapause development days DPSW Summer or winter diapause ('0' for winter diapause, '1'for summer diapause) TTCS Cold stress temperature threshold THCS Cold stress accumulation rate TTHS Heat stress temperature threshold THHS Heat stress accumulation rate SMDS Dry stress soil moisture threshold HDS Dry stress accumulation rate SMWS Wet stress soil moisture threshold Estimating the size of land areas of suitable climate projected by the

Fig 1

 1 Fig 1 Climatic suitability (EI) for the locust in the current and future climate scenarios projected using CLIMEX™. a)Qinghai-Tibet Plateau, current climate; b)Qinghai-Tibet Plateau, 2020; c)Qinghai-Tibet Plateau, 2050; d)Qinghai-Tibet Plateau, 2100. □, unsuitable (0.00); , marginal (0.01-10); , suitable (10-20); , optimal (20.00+).

Acknowledgements

The authors thank Prof. Stauffer at Pennsylvania State University for his valuable suggestions and modifying the draft. We would like to thank Yan Fan, Zhimei Li and all other members of the Plant Quarantine Laboratory of China Agricultural University (CAUPQL); Also we would like to thank Dezhi Li and other members of Key Laboratory for Biological Control of Ministry of Agriculture, Institute of Plant Protection, Chinese Academy of Agricultural Sciences.

This study received financial supports from the Special Fund for Agro-scientific Research in the Public Interest (201003079).