CN103761447B - Planting area remote sensing confirming method for grain subsidy accounting of prefecture level and county level - Google Patents

Abstract

The invention discloses a planting area remote sensing confirming method for grain subsidy accounting of the prefecture level and the county level. The method comprises the following steps: P1, conducting sample investigation on the planting area of grain corps in a province; P2, selecting and preprocessing a remote-sensing image; P3, conducting remote sensing recognition on the grain crops; P4, processing mixed pixels; P5, establishing a grain crop planting area estimation model, and estimating the planting area of large quantity of the grain crops in district-and-county-level administrative units.

Description

The cultivated area remote sensing adjusted for districts and cities and county's two-stage grain subsidy determines method

Technical field

The present invention relates to remote-sensing monitoring method, more particularly to a kind of grain adjusted for districts and cities and county's two-stage grain subsidy Cultivated area remote sensing estimation method.

Background technology

Grain direct payment policy is that mostly important government of China supports agriculture one of policy in recent years, it is intended to " further plus Power Huinong, policy of supporting agriculture dynamics, preferably transfer peasant and plant grain enthusiasm, protection grain producer interests ".With to peasant Subsidy fund is continuously increased, and government is also required to strengthen to the management for subsidizing fund.Each province of current China is started with grain The subsidy policy on cultivated area basis.In order to implement the subsidy of Farming, it is necessary to districts and cities and the grain of administration cell at county level Cultivated area is adjusted.Statistical report mode is presently mainly used, foodstuff planting area is reported by each county, then by provincial The finance department adjusts the subsidy amount of money according to the foodstuff planting area for reporting.

Estimate that crop acreage is the Main Means of current statistical department by statistical, but due to sample cloth If, ground investigation and statistical error, and subjective aspect factor, the quality of data is difficult to ensure that.And due to the data be for The grain subsidy amount of money in each county oneself is calculated, therefore there is obvious false phenomenon, state revenue and expenditure fund is there is larger wave Take phenomenon.Therefore a kind of objective, reliable monitoring method must be developed, and remote sensing technology is because it is on time and space yardstick Superiority, possesses advantage in terms of large plant of grain crops area is quick and precisely monitored.Typically exist yet with grain subsidy Carried out on the level administration cell of district, thus require to have in area's County Scale with accuracy (being better than 85%) higher, make profit Monitored with remote sensing technique and large plant of grain crops area of district level is accurate in big region to turn into a difficulty larger and anxious Need to solving for task.

The content of the invention

For the deficiencies in the prior art, districts and cities and county's two-stage grain to be used for it is an object of the invention to provide one kind Subsidize the cultivated area remote sensing adjusted and determine method, can effectively reduce and obtain crop-planting face by ground investigation interior on a large scale The cost of product information, while the subjective impact produced by artificial participation during being also effectively prevented from ground investigation, is The execution of China's food subsidy policy provides objective, accurate Back ground Information.

The technical proposal of the invention is realized in this way：A kind of growing surface adjusted for districts and cities and county's two-stage grain subsidy Product remote sensing determines method, comprises the following steps：

Step P1：Cultivated area to cereal crops in target province domain is sampled investigation；

Step P2：Remote sensing image is selected and pretreatment；

Step P3：Cereal crops remote sensing recognition；

Step P4：The treatment of mixed pixel；

Step P5：The crop area estimation model of cereal crops is set up, the grain to districts and cities and county's two-stage administration cell is made The cultivated area of thing is estimated.

The above-mentioned cultivated area remote sensing adjusted for districts and cities and county's two-stage grain subsidy determines method, in step P1：With Based on research area's history crop acreage and planting proportion structure, the area sampling frame of province domain scope, and root are set up According to the pattern of farming of cereal crops, expectation quality and Calculation of Reliability sample size, and cloth specimen is carried out using system cloth sample loading mode, had Body comprises the following steps：

Step P101：The rectangle grid of covering province domain scope, rectangle Grid size are set up using GIS-Geographic Information System instrument Size is：Length and width is 5 kilometers, deletes the rectangle grid unit across provincial boundaries, according to from top to bottom, from the left and right side Order is numbered to each rectangle grid unit so that each rectangle grid unit has unique numbering；

Step P102：The cultivated area and planting proportion structure of the cereal crops of research on utilization area recent five years, calculate sample Amount, Sample Size Calculation is as follows：

N=n₀+n₀* 5%；

In formula：n₀It is initial sample size, n is amendment sample size；T is sampling probability, when confidence level is 95%, t etc. In 1.96；R is relative error, and r=5% is taken in research；It is population mean, that is, studies the kind of the cereal crops of area's recent five years Plant area average；S²It is population variance, that is, studies the cultivated area variance of the cereal crops of area's recent five years；

Step P103：Using system cloth quadrat method, ground investigation sample prescription is selected, specifically include following steps：(1) with total Rectangular grid netting index obtains cloth specimen interval divided by sample size；(2) randomly select one between 0 and cloth specimen interval numerical value between with Machine number, as the initial value of cloth specimen, that is, extracts first sample；(3) the rectangle grid numbering according to a upper sample, cumulative one Individual cloth specimen interval, obtains a rectangle grid unit numbering for new sample；(4) repeat step (3), until terminating；

Step P104：The proportion of crop planting classification on every piece of arable land in manual research each sample rectangle grid；

Step P105：The crop proportion of the every kind of crop of statistical inference, and according to cultivated area statistical estimation the whole province of the whole province Each crop acreage of the whole province of cultivated area, the i.e. methods of sampling estimation；Infer that formula is as follows：

Total cultivated area=the proportion of crop planting of crops is into several * total area under cultivations.

The above-mentioned cultivated area remote sensing adjusted for districts and cities and county's two-stage grain subsidy determines method, in step P2：Root According to the crop phenological period, it is determined that the environment star ccd data of suitable phase is used as data source, and selected environment star ccd data is entered Row pretreatment, specifically includes following steps：

Step P201：With reference to research area's history phenology data, respectively 30 days after crop sowing time, the jointing stage, florescence 4 scape images are obtained with preharvest period；

Step P202：With 1:100000 topographic maps are to refer to image, and environment star CCD images are entered using quadratic polynomial method Row geometric accurate correction, control errors are within 1 pixel；

Step P203：Radiation calibration treatment is carried out to environment star CCD images, calibration is carried out using equation below：

In formula, DN is the gray value of remote sensing image, A and L₀For the environment disaster reduction constellation that Chinese Resources satellite hub is announced The in-orbit Absolute Radiometric Calibration Coefficients of each load of A/B stars, d is day ground unit distance, and ESUN is air top layer solar irradiance, θ is solar zenith angle；

Step P204：It is the ratio vegetation index RVI of calculating optical remote sensing image, normalized differential vegetation index NDVI, enhanced Vegetation index EVI：

Wherein, ρ_BlueIt is environment star ccd data indigo plant wave band reflectivity, ρ_NirIt is near infrared reflectivity, ρ_RedFor red wave band is anti- Penetrate rate.L is adjusting parameter, C₁、C₂It is Atmospheric corrections parameter, they can reduce the effect of background and air, for environment star For image, L, C₁、C₂Empirical value be respectively 1,6 and 7.5；

Step P205：By NDVI data together with 4 band overlappings of environment star ccd data, generation has 5 wave bands Composograph；

Step P206：Cultivation area is extracted：The arable land distributed data of the whole province is extracted using the whole province's arable land database, with complete It is mask to save arable land distributed data, and the environment star chart picture to synthesizing does mask process, arable land cut out, as the base of subsequent treatment Plinth data.

The above-mentioned cultivated area remote sensing adjusted for districts and cities and county's two-stage grain subsidy determines method, in step P3：Profit With maximum likelihood classification or support vector machine classification method, with step P104 ground sample investigation results as training sample, to ring Border star CCD images are classified, and recognize different cereal crops, specifically include following steps：

Step P301：Ground investigation sample vector data are called in, two faces are divided into using GIS-Geographic Information System instrument The equal part of product, will sample be divided into two Sample Storehouses；Different cereal crops polygons is extracted respectively, and to every Individual polygon is numbered so that each polygon has unique identifying number；According to polygonal unique identifying number, grain is made The polygon of thing is divided into 2 parts, and odd number record part 1, even number is recorded as part 2, while deleting area is less than 8100 Square metre polygon；

Step P302：Any part in two Sample Storehouses obtained by selecting step P301, as crops remote sensing point The training sample database of class, therefrom chooses that plot is larger and polygon of color stable, as spectrum training set, then using spectrum Training tool, trains the spectral signature of different crops；

Step P303：Image classification is carried out using maximum likelihood method or SVMs, each pixel is assigned to certain agriculture Crop, finally gives the Crop classification result of image, and classification results are recoded, a kind of corresponding pixel of crop Value is set to 1, and other crops are set to 0, obtains the pure pixel distributed data collection of the crop；This step is repeated to every kind of crop, respectively To the pure pixel distributed data collection of every kind of crop, so as to obtain the distribution map of different crops.

The above-mentioned cultivated area remote sensing adjusted for districts and cities and county's two-stage grain subsidy determines method, in step P4： On the basis of step P3 classification, for the mixed pixel of multiple kinds of crops interlaced area, the decomposition mould based on spectral signature is set up Type, calculates the area percentage of various crops in different pixels, specifically includes following steps：

Step P401：Image is obtained to P206 carries out unsupervised classification, obtains the planting area of different crops, and removes The pure picture dot region that step P303 is obtained, obtains mixed pixel distributed area；

Step P402：The spectral signature (average value) of Different Crop is read from the spectrum training characteristics of step P3, is calculated Every kind of crops and bare area RVI, NDVI, the average of EVI, set up linear spectral unmixing model, and computing formula is：

Wherein, n is proportion of crop planting quantity, p_iA kind of pixel scale shared by crops, p_lPicture shared by the bare area First ratio, RVI_i、NDVI_i、EVI_iIt is average RVI values, NDVI values, the EVI values, RVI of the training sample of the crops_l、NDVI_l、 EVI_lThe average RVI values of the training sample of bare area, NDVI values, EVI values, RVI_p、NDVI_p、EVI_pFor correspondence pixel RVI values, NDVI values, EVI values；

Step P403：Using linear spectral unmixing model, the Different Crop cultivated area ratio of different pixels is solved.

The above-mentioned cultivated area remote sensing adjusted for districts and cities and county's two-stage grain subsidy determines method, in step P5：With Based on step P1, step P3 and step P4, plant of grain crops area appraising model is set up, to the grain of district level administration cell Food crop cultivated area estimated, specifically includes following steps：

Step P501：According to the area ratio of different pixels, different grains in different districts and cities and county's two-stage administration cell are counted The pel data of food crop, and the hair plantation of various crops in different districts and cities and county's two-stage administration cell is calculated according to Pixel size Area；

Step P502：Using the Cultivation land coefficient in different districts and cities and county's two-stage administration cell, (arable land accounts for the ratio of figure spot area Example), estimate the net cultivated area of various crops in each districts and cities and county's two-stage administration cell：

Net cultivated area=hair cultivated area * Cultivation land coefficients；

Step P503：Use regression estimation model, the kind of the crop that the crop proportion of the crop that sampling is obtained and classification are obtained Number mutually amendment is planted into, the cultivated area of prefectures and cities and the county various cereal crops of two-stage administration cell is calculated；

Step P504：Statistical summaries obtain the cultivated area of the various cereal crops of the whole province.

The invention has the advantages that：

1. ability of the remote sensing technology in terms of crop identification on a large scale is made full use of, it is interior on a large scale so as to effectively reduce The cost of crops planting area information is obtained by ground investigation, while because of artificial participation during also effectively prevent ground investigation And the subjective impact for producing；

2. as data source, combined ground sample investigation data pass through No. 1 satellite data of environment with China's independent research Each wave band reflectivity of ccd data and normalized differential vegetation index, analyze spectrum, the Texture eigenvalue of large cereal crops, and utilize The methods such as maximum likelihood classification, SVMs recognize large cereal crops, on this basis, with reference to second tune of land resources Arable land database and ridge coefficient data storehouse are looked into, large plant of grain crops area appraising model is set up, and rower is entered to model Fixed, by provincial administration cell refine to reporting unit at county level by the crops planting area monitoring method based on remote sensing for developing Administration cell, the large plant of grain crops area filled up in district level reporting unit estimates the blank in field, is that grain is mended The basic crop area estimation of patch policy provides important technical support, so that for the execution of China's food subsidy policy is provided Objective, accurate Back ground Information；

3. No. 1 satellite of environment has room and time resolution ratio higher, and MODIS numbers are substantially better than in spatial resolution According to, Landsat TM data are substantially better than in temporal resolution, while its breadth is wider, ensure that crops planting area industry The need for businessization is monitored.

4. the crop area estimation of the cereal crops such as wheat, corn and soybean, paddy rice, potato is applied to.

Brief description of the drawings

Fig. 1 early rice sample investigation areas and the classification area regression analysis of environment star；

Crops planting area remote sensing estimation method flow chart of Fig. 2 present invention for grain subsidy；

Fig. 3 semilate rice sample investigation areas and the classification area regression analysis of environment star；

Fig. 4 late rice sample investigation areas and the classification area regression analysis of environment star.

Specific embodiment

With reference to accompanying drawing, the present invention is described further：

As shown in Fig. 2 the present invention determines method bag for the cultivated area remote sensing that districts and cities and county's two-stage grain subsidy are adjusted Include following steps：

Step P1：Cultivated area to cereal crops in target province domain is sampled investigation.

Based on studying area's history crop acreage and planting proportion structure, the area sampling of province domain scope is set up Framework, and pattern of farming according to cereal crops, expectation quality and Calculation of Reliability sample size, and entered using system cloth sample loading mode Row cloth specimen, specifically includes following steps：

Step P101：Rectangle grid is produced using the GIS-Geographic Information System such as ArcGIS (GIS) instrument, rectangle Grid size is big It is small to be：Length and width is 5 kilometers, it is desirable to covering province domain scope.The rectangle grid unit across provincial boundaries is deleted, then to every According to from top to bottom, from a left side, the order on the right side is numbered individual rectangle grid unit so that each rectangle grid unit has only One numbering.

Step P102：The cultivated area and planting proportion structured data of the cereal crops of research on utilization area recent five years, calculate Sample size, Sample Size Calculation is as follows：

N=n₀+n₀* 5%；

In formula：n₀It is initial sample size, n is amendment sample size；T is sampling probability, when confidence level is 95%, t etc. In 1.96；R is relative error, and 5% is taken in research；It is population mean, i.e., the cultivated area of the cereal crops of 5 years is average；S²It is Population variance, that is, study the cultivated area variance of the cereal crops of area's recent five years.In order to ensure sampling precision, sample size is in reality 5% is raised in the implementation procedure of border.

Step P103：Using system cloth quadrat method, ground investigation sample prescription is selected, concrete operations are：(1) with total rectangular grid Netting index obtains cloth specimen interval divided by sample size；(2) random number between 0 and cloth specimen interval numerical value is randomly selected, is made It is the initial value of cloth specimen, that is, extracts first sample；(3) the rectangle grid numbering according to a upper sample, add up a cloth specimen Interval, obtains a rectangle grid unit numbering for new sample；(4) repeat step (3), until terminating.

Step P104：The crop-planting classification on every piece of arable land in manual research each sample rectangle grid.Specific behaviour Work is：Using the figure spot figure in GIS-Geographic Information System (GIS) treatment sample grid, sketched out each plot with differential GPS measure Border, corresponding numbering is then filled according to intramassif Crop Group, specific coding is as follows：

1- winter wheat；2- mono- season rice/semilate rice；3- corns；4- spring wheat；5- early rice；6- late rices；7- soybean；9- sorghums； 10- potatos；11- other beans；12- lies fallow；13- is no longer cultivated.

Step P105：The crop proportion of the every kind of crop of statistical inference, and according to cultivated area statistical estimation the whole province of the whole province Each crop acreage of the whole province of cultivated area, the i.e. methods of sampling estimation；Infer that formula is as follows：

Total cultivated area=the proportion of crop planting of crops is into several * total area under cultivations.

Step P2：Remote sensing image is selected and pretreatment.

According to the crop phenological period, it is determined that the environment star ccd data of suitable phase is used as data source, and the environment star to selecting Ccd data is pre-processed, and is specifically included：

Step P201：With reference to research area's history phenology data, respectively 30 days after crop sowing time, the jointing stage, florescence 4 scape images are obtained with preharvest period；

The crop phenological period for estimating area (monitoring section or province) according to target selects suitable environment star CCD images, natural seeding 1 month image to before gathering in can be used after kind, in order to improve nicety of grading, it is desirable to which at least 2 phase images, interval time should surpass Spend 20 days.

Step P202：With 1:100000 topographic maps are to refer to image, and geometric accurate correction is carried out to environment star CCD images, correction Method uses quadratic polynomial method, and method for resampling uses closest Furthest Neighbor；Control errors：Plains region requires error not More than 1 pixel, mountain area error is no more than 1.5 pixels.

Step P203：Radiation calibration treatment is carried out to environment star CCD images, calibration is carried out using equation below：

In formula, DN is the gray value of remote sensing image, A and L₀For the environment disaster reduction constellation that Chinese Resources satellite hub is announced The in-orbit Absolute Radiometric Calibration Coefficients of each load of A/B stars, d is day ground unit distance, and ESUN is air top layer solar irradiance, θ is solar zenith angle；

Step P204：Ratio calculated vegetation index RVI, normalized differential vegetation index NDVI, enhancement mode meta file EVI, meter Calculate formula as follows：

Wherein, ρ_BlueIt is environment star ccd data indigo plant wave band reflectivity, ρ_NirIt is near infrared reflectivity, ρ_RedFor red wave band is anti- Penetrate rate.L is adjusting parameter, C₁、C₂It is Atmospheric corrections parameter, they can reduce the effect of background and air, for environment star For image, L, C₁、C₂Empirical value be respectively 1,6 and 7.5；

Step P205：By 4 band overlappings of normalized differential vegetation index (NDVI) data and environment star ccd data one Rise, generation has 5 composographs of wave band.

Step P206：Cultivation area is extracted：Ploughing for the whole province is extracted using target estimation area (such as the whole province) arable land database Ground distributed data, distributed data is ploughed as mask with target estimation area (such as the whole province), and the environment star chart picture to synthesizing is done at mask Reason, cuts out arable land, as the basic data of subsequent treatment.

Step P3：Cereal crops remote sensing recognition.

Using maximum likelihood classification or support vector machine classification method, with ground sample investigation as training sample, to environment Star CCD images are classified, and recognize different cereal crops, specifically include following steps：

Step P301：Ground investigation sample vector data are called in, two faces are divided into using GIS-Geographic Information System instrument The equal part (being divided into two Sample Storehouses) of product, respectively extracts the polygon of different cereal crops, and many to each Side shape is numbered so that each polygon has unique identifying number.

According to polygonal unique identifying number, the polygon of cereal crops is divided into 2 parts (such as：Odd number is recorded as the 1st Part, even number is recorded as part 2)；Area is deleted simultaneously less than 8100 square metres of polygons of (9 environment star image units)【Face Product is too small without representativeness, cannot use】.

Step P302：Any portion in two Sample Storehouses obtained by selecting step P301, as crops Classification in Remote Sensing Image Training sample database, therefrom choose that plot is larger and polygon of corresponding environment sing data color stable, trained as spectrum Collection, then using spectrum training tool, calculates maximum, minimum value, average, the variance of each wave band of polygon correspondence pixel Deng as the spectral signature of different crops.

Step P303：Image classification is carried out using maximum likelihood method or SVMs, by each as being assigned to certain farming Thing.

Classification results are recoded, a kind of corresponding pixel value of crop is set to 1, other crops are set to 0, are somebody's turn to do The pure pixel distribution of crop.This step is repeated to every kind of crop, the pure pixel distribution of every kind of crop is respectively obtained, so as to obtain different agricultures The distribution map of crop.

Step P4：The treatment of mixed pixel.

Decomposition of Mixed Pixels is carried out to multiple kinds of crops interlaced area, the linear unmixed model based on spectral signature is set up, The area percentage of various crops in different pixels is calculated, following steps are specifically included：

Step P401：Environment sing data to synthesizing carries out unsupervised classification, post-classification comparison, pixel cluster group difference Various ground class is divided into, is represented with corresponding code.Non-supervised classification obtains proportion of crop planting region, removes supervised classification The pure picture dot region that method is obtained, obtains mixed pixel distributed area.

Step P402：According to the spectral signature in step P3, calculate every kind of crops and bare area RVI, NDVI, EVI it is equal Value, using linear spectral unmixing model, sets up equation, the ratio in calculating pixel shared by the crops, and computing formula is：

Wherein, n is proportion of crop planting quantity, p_iA kind of pixel scale shared by crops, p_lPicture shared by the bare area First ratio, RVI_i、NDVI_i、EVI_iIt is average RVI values, NDVI values, the EVI values, RVI of the training sample of the crops_l、NDVI_l、 EVI_lThe average RVI values of the training sample of bare area, NDVI values, EVI values, RVI_p、NDVI_p、EVI_pFor correspondence pixel RVI values, NDVI values, EVI values；

Step P403：Using linear spectral unmixing model, the Different Crop cultivated area ratio of different pixels is solved.

Step P5：Based on step P1, step P3 and step P4, plant of grain crops area appraising model is set up, it is right The plant of grain crops area of districts and cities and county's two-stage administration cell estimated, specifically includes following steps：

Step P501：Similar atural object data to being obtained in step P3 and step P4 are inlayed, overlapping region pixel value Selection maximum, obtains the crop distribution comprising pure pixel and mixed pixel.According to the area ratio of different pixels, statistics is different The pel data of different plant of grain crops in districts and cities and county's two-stage administration cell, and according to Pixel size calculate different districts and cities and The hair cultivated area of various crops in county's two-stage administration cell.

Step P502：Using the Cultivation land coefficient in prefectures and cities and county's two-stage administration cell, (arable land accounts for the ratio of figure spot area Example), the net cultivated area of various crops in estimation districts and cities and county's two-stage administration cell, i.e.,：

Net cultivated area=hair cultivated area * Cultivation land coefficients,

Step P503：By prefectures and cities and county's two-stage administration cell crop area summation, the whole province is obtained using remote sensing technique Different crops cultivated area, and the whole province's cultivated area is combined, calculate different crops crop proportion across the entire province

The crop area that classification is obtained in each sample prescription is counted, proportion of crop planting into number is calculated in the range of sample prescriptionThen The crop area that the crop area obtained to ground investigation in each sample prescription is obtained with classification carries out regression analysis, calculates back Return equation, obtain regression coefficient b.

Such as in Fig. 1, b=1.105；Such as in figure 3, b=0.856；Such as in fig. 4, b=0.917.In Fig. 1,3,4, test The size for demonstrate,proving sample prescription is 1*1KM, therefore according to this result, we can draw our evaluation method 1*1KM's Yardstick has had precision very high, for thicker County Scale, region yardstick, excellent layout for, estimation precision is obvious Can be what is fully ensured that.

Field is estimated in remote sensing recognition, because southern area plot relatively crushes, pattern of farming is complicated, and Rice Cropping point It it is early, middle and late three season, semilate rice and early late rice have the situation that growth time overlaps, and its plantation situation is substantially multiple compared with other crops It is miscellaneous, therefore, Monitoring of Paddy Rice Plant Area is that being most difficult to of generally acknowledging in cereal crops is identified estimation, by Fig. 1, Fig. 3 and Fig. 4 institute Show, this case technology scheme can be very accurately estimated the cultivated area of early rice, semilate rice and late rice, thus this example Technical scheme can also carry out more accurately estimation to other plant of grain crops area in addition to paddy rice.

Each crop acreage of the whole province estimated using the methods of sampling and sorting technique, sets up different crops growing surface Product regression estimation model：

Wherein,The crop-planting of regression estimation Modifying model into number,The crop-planting obtained for ground investigation into Number, b is correction factor,Be sorting technique statistics proportion of crop planting across the entire province into number,It is sorting technique statistics Proportion of crop planting is into number in the range of sample prescription.

Step P504：The crop-planting corrected is into numberAfterwards, by prefectures and cities and county's two-stage administration cell cultivated area Data calculate prefectures and cities and county's two-stage administration cell different crops cultivated area.

Above-described embodiment is only intended to clearly illustrate the invention example, and not has to the invention The restriction of body implementation method.For those of ordinary skill in the field, can also make on the basis of the above description The change or variation of other multi-forms.There is no need and unable to be exhaustive to all of implementation method.It is all of the invention Spirit and principle within extend out it is any obvious change or variation still in the invention claim guarantor Among shield scope.

Claims (2)

1. the cultivated area remote sensing for being used for districts and cities and county's two-stage grain subsidy accounting determines method, it is characterised in that including as follows Step：

Step P1：Cultivated area to cereal crops in target province domain is sampled investigation；

Step P2：Remote sensing image is selected and pretreatment；

Step P3：Cereal crops remote sensing recognition；

Step P4：The treatment of mixed pixel；

Step P5：The crop area estimation model of cereal crops is set up, to the cereal crops of districts and cities and county's two-stage administration cell Cultivated area is estimated；

In step P1：Based on monitoring province crop acreage and planting proportion structure history data, province domain is set up The area sampling frame of scope, and pattern of farming according to cereal crops, expectation quality and Calculation of Reliability sample size, and use System cloth sample loading mode carries out cloth specimen, specifically includes following steps：

Step P101：The rectangle grid of covering province domain scope, rectangle grid size are set up using GIS-Geographic Information System instrument For：Length and width is 5 kilometers, deletes the rectangle grid unit across provincial boundaries, according to from top to bottom, from the left and order on the right side Each rectangle grid unit is numbered so that each rectangle grid unit has unique numbering；

Step P102：The cultivated area and planting proportion structured data of the cereal crops of research on utilization area recent five years, calculate sample Amount, Sample Size Calculation is as follows：

$n_0={\left(\frac{t}{r}\right)}^2\frac{S^2}{{\stackrel{\‾}{Y}}^2},$

N=n₀+n₀* 5%；

In formula：n₀It is initial sample size, n is amendment sample size；T is sampling probability, and when confidence level is 95%, t is equal to 1.96；R is relative error, and r=5% is taken in research；It is population mean, that is, studies the growing surface of the cereal crops of area's recent five years Product is average；S²It is population variance, that is, studies the cultivated area variance of the cereal crops of area's recent five years；

Step P103：Using system cloth quadrat method, ground investigation sample prescription is selected, specifically include following steps：(1) with total rectangle Grid number obtains cloth specimen interval divided by sample size；(2) random number between 0 and cloth specimen interval numerical value is randomly selected, As the initial value of cloth specimen, that is, extract first sample；(3) the rectangle grid numbering according to a upper sample, add up a cloth Sample is spaced, and obtains a rectangle grid unit numbering for new sample；(4) repeat step (3), until terminating；

Step P104：The proportion of crop planting classification on every piece of arable land in each sample rectangle grid is investigated by special troop；

Step P105：The crop proportion of the every kind of crop of statistical inference, and according to the plantation of cultivated area statistical estimation the whole province of the whole province Area, infers that formula is as follows：

Total cultivated area=the proportion of crop planting of crops is into several * total area under cultivations；

In step P2：According to the crop phenological period, it is determined that the environment star ccd data of suitable phase is used as data source, and to selected Environment star ccd data pre-processed, specifically include following steps：

Step P201：With reference to research area's history phenology data, respectively 30 days, jointing stage, florescence and receipts after crop sowing time Obtain early stage and obtain 4 scape images；

Step P202：With 1:100000 topographic maps are to refer to image, environment star CCD images are carried out using quadratic polynomial method several What fine correction, control errors are within 1 pixel；

Step P203：Radiation calibration treatment is carried out to environment star CCD images, calibration is carried out using equation below：

$\mathrm{\ρ}=(\frac{DN}{A}+L_0)*\frac{\mathrm{\π}*d^2}{ESUN*cos\mathrm{\θ}},$

In formula, DN is the gray value of remote sensing image, A and L₀For the environment disaster reduction constellation A/B stars that Chinese Resources satellite hub is announced Each in-orbit Absolute Radiometric Calibration Coefficients of load, d is day ground unit distance, and ESUN is air top layer solar irradiance, and θ is for too Positive zenith angle；

Step P204：The ratio vegetation index RVI of calculating optical remote sensing image, normalized differential vegetation index NDVI, enhanced vegetation Index E VI：

$RVI=\frac{{\mathrm{\ρ}}_{Nir}}{{\mathrm{\ρ}}_{\mathrm{Re}d}};$

$NDVI=\frac{{\mathrm{\ρ}}_{Nir}-{\mathrm{\ρ}}_{\mathrm{Re}d}}{{\mathrm{\ρ}}_{Nir}+{\mathrm{\ρ}}_{\mathrm{Re}d}};$

$EVI=(1.5+L)\×\frac{{\mathrm{\ρ}}_{Nir}-{\mathrm{\ρ}}_{\mathrm{Re}d}}{L+{\mathrm{\ρ}}_{Nir}+C_1\×{\mathrm{\ρ}}_{\mathrm{Re}d}+C_2\×{\mathrm{\ρ}}_{Blue}}$

Wherein, ρ_BlueIt is environment star ccd data indigo plant wave band reflectivity, ρ_NirIt is near infrared reflectivity, ρ_RedIt is red wave band reflectivity, L It is adjusting parameter, C₁、C₂It is Atmospheric corrections parameter, they can reduce the effect of background and air, comes for environment star image Say, L, C₁、C₂Empirical value be respectively 1,6 and 7.5；

Step P205：By NDVI data together with 4 band overlappings of environment star ccd data, generation has 5 conjunctions of wave band Into image；

Step P206：Cultivation area is extracted：The arable land distributed data of the whole province is extracted using the whole province's arable land database, is ploughed with the whole province Ground distributed data is mask, and the environment star chart picture to synthesizing does mask process, arable land cut out, as the basic number of subsequent treatment According to；

In step P3：Using maximum likelihood classification or support vector machine classification method, investigated with step P104 ground samples and tied Fruit is training sample, and environment star CCD images are classified, and recognizes different cereal crops, specifically includes following steps：

Step P301：Ground investigation sample vector data are called in, two area phases are divided into using GIS-Geographic Information System instrument Deng part, will sample be divided into two Sample Storehouses；Different cereal crops polygons is extracted respectively, and it is many to each Side shape is numbered so that each polygon has unique identifying number；According to polygonal unique identifying number, by cereal crops Polygon is divided into 2 parts, and odd number is recorded as part 1, and even number is recorded as part 2, while it is flat less than 8100 to delete area The polygon of square rice；

Step P302：Any part in two Sample Storehouses obtained by selecting step P301, as crops Classification in Remote Sensing Image Training sample database, therefrom chooses plot greatly and the polygon of color stable, as spectrum training set, then trains work using spectrum Tool, trains the spectral signature of different crops；

Step P303：Image classification is carried out using maximum likelihood method or SVMs, each pixel is assigned to certain crops, The Crop classification result of image is finally given, and classification results are recoded, a kind of corresponding pixel value of crop is set It is 1, other crops are set to 0, obtains the pure pixel distributed data collection of the crop；This step is repeated to every kind of crop, is respectively obtained every The pure pixel distributed data collection of crop is planted, so as to obtain the distribution map of different crops；

In step P4：On the basis of step P3 classification, for the mixed pixel of multiple kinds of crops interlaced area, foundation is based on The decomposition model of spectral signature, calculates the area percentage of various crops in different pixels,

Specifically include following steps：

Step P401：Image is obtained to P206 carries out unsupervised classification, obtains the planting area of different crops, and removal step The pure picture dot region that P303 is obtained, obtains mixed pixel distributed area；

Step P402：The spectral signature average value of Different Crop is read from the spectrum training characteristics of step P3, every kind of agriculture is calculated Crop and bare area RVI, NDVI, the average of EVI, set up linear spectral unmixing model, and computing formula is：

$\underset{i=1}{\overset{n}{\Σ}}{p}_i+p_l=1$

$\underset{i=1}{\overset{n}{\Σ}}{p}_i*{\mathrm{RVI}}_i+p_l*{\mathrm{RVI}}_l={\mathrm{RVI}}_p$

$\underset{i=1}{\overset{n}{\Σ}}{p}_i*{\mathrm{NDVI}}_i+p_l*{\mathrm{NDVI}}_l={\mathrm{NDVI}}_p$

$\underset{i=1}{\overset{n}{\Σ}}{p}_i*{\mathrm{EVI}}_i+p_l*{\mathrm{EVI}}_l={\mathrm{EVI}}_p$

Wherein, n is proportion of crop planting quantity, p_iA kind of pixel scale shared by crops, p_lPixel ratio shared by the bare area Example, RVI_i、NDVI_i、EVI_iIt is average RVI values, NDVI values, the EVI values, RVI of the training sample of the crops_l、NDVI_l、EVI_l The average RVI values of the training sample of bare area, NDVI values, EVI values, RVI_p、NDVI_p、EVI_pIt is RVI values, the NDVI of correspondence pixel Value, EVI values；

Step P403：Using linear spectral unmixing model, the Different Crop cultivated area ratio of different pixels is solved.

2. the cultivated area remote sensing adjusted for districts and cities and county's two-stage grain subsidy according to claim 1 determines method, Characterized in that, in step P5：Based on step P1, step P3 and step P4, plant of grain crops area estimation is set up Model, the plant of grain crops area to districts and cities and county's two-stage administration cell estimates, specifically includes following steps：

Step P501：According to the area ratio of different pixels, count different grains in different districts and cities and county's two-stage administration cell and make The pel data of thing, and the hair growing surface of various crops in different districts and cities and county's two-stage administration cell is calculated according to Pixel size Product；

Step P502：Using the Cultivation land coefficient in different districts and cities and district, estimate each in each districts and cities and county's two-stage administration cell Plant the net cultivated area of crop：

Net cultivated area=hair cultivated area * Cultivation land coefficients；

Step P503：Use regression estimation model, the plantation of the crop that the crop proportion of the crop that sampling is obtained is obtained with classification into Number mutually amendment, calculates the cultivated area of various cereal crops in prefectures and cities and county's two-stage administration cell；

Step P504：Statistical summaries obtain the cultivated area of the various cereal crops of the whole province.