CN105046188A - MODIS mixed pixels decomposition forest information extraction method - Google Patents

Abstract

The invention provides an MODIS mixed pixels decomposition forest information extraction method, being integrated with the standard images of an MOD09A1 surface reflectance product and the standard images of an MOD13Q1 vegetation index product. The MODIS mixed pixels decomposition forest information extraction method is characterized in that the images of different MODIS products are provided with different forest information; and however, as the resolution ratio is low and large amount of mixed pixels exist in the images, the forest information is extracted through an improved end-member purification method and adoption of a linear mixed pixel decomposition model. The MODIS mixed pixels decomposition forest information extraction method of the invention has the advantages of being able to quickly extract the forest cover information in wide range, improving the forest mapping precision and the forest type identification precision, wherein only the MOD09A1 surface reflectance product and the MOD13Q1 vegetation index product are used during the information extraction process, thus obtaining ideal effect and being simple and practicable.

Description

MODIS Decomposition of Mixed Pixels forest information extracting method

Art

The present invention relates to the technology of a kind of MODIS remote sensing image forest information extraction, Decomposition of Mixed Pixels especially can be utilized to improve forest cover drawing and Forest Types accuracy of identification sensor information extracting method.

Background technology

The remote sensing image of current middle high-resolution, although measuring accuracy is higher, the price comparison of data is expensive, acquisition is not easy; Low resolution remote sensing image coverage is wider; and expense is cheaper; as MODIS image; because its most high spatial resolution is 250m, and MODIS remote sensing image has very high spectral radiance precision, can Free Acquisition; but there is a large amount of mixed pixel; when the classification of atural object, usually can cause certain error, utilize classic method to be difficult to obtain desirable nicety of grading.

Summary of the invention

Need of production and the poor deficiency of Forest Types recognition capability can not be met to overcome existing MODIS image forest information extracting method nicety of grading, the invention provides a kind of forest information extracting method, the feature that the method can effectively utilize MODIS data spectrum Electrodynamic radiation splendid, can improve again the deficiency of MODIS spatial resolution, and the MOD09A1 Reflectivity for Growing Season product that MODIS image can be utilized quickly and easily to provide and MOD13Q1 vegetation index product effectively improve forest cover drawing and Forest Types accuracy of identification.

In order to realize above-mentioned technical purpose, technical scheme of the present invention is, a kind of MODIS Decomposition of Mixed Pixels forest information extracting method, comprises the steps:

Step one: data acquisition step, obtains forest remote sensing image;

Step 2: pre-treatment step, carries out pre-service to the remote sensing image that step one obtains, obtains pretreated remote sensing image data;

Step 3: Forest Types phenology difference characteristic analytical procedure, finds out the vegetation index period obviously can distinguishing each Forest Types;

Step 4: end-member composition purification step, obtains terminal pixel component;

Step 5: end-member composition spectral reflectivity asks for step, obtains each reflectance;

Step 6: linear spectral unmixing step, asks for each atural object distribution plan;

Step 7: forest information extracting step, obtains each atural object distribution plan.

Described method, in described step 2, the pre-service of remotely-sensed data comprises: carry out projection transform to remotely-sensed data, atmospheric correction, image mosaic, eliminates black stripe, geometry correction, and study area is extracted.

Described method, the step that described step 3 comprises is:

GPS is utilized to obtain the coordinate information of typical sampling point, GIS software is utilized to obtain the gray-scale value of typical sampling point on remote sensing image, build vegetation index time series section curve, the data that the vegetation index remotely-sensed data obviously can distinguishing each Forest Types then found out by phenology analysis and MOD09A1 Reflectivity for Growing Season product combined together are obtained.

Described method, described typical sampling point comprises coniferous forest, broad-leaf forest, bamboo grove, spinney, waters, arable land and construction land.

Described method, the data that the described vegetation index remotely-sensed data obviously can distinguishing each Forest Types found out by phenology analysis and MOD09A1 Reflectivity for Growing Season product combined together are obtained comprise the following steps:

First end member purification model is built, utilize decision tree classification computer graphics sorting technique, according to the feature of each Forest Types phenology difference, utilize the vegetation index remotely-sensed data in period obviously can distinguishing each Forest Types found out, and utilize the vegetation index gray-scale value can distinguishing each atural object as the threshold value of Decision-Tree Classifier Model: namely

When the threshold value of the 8th phase NDVI be less than 0.6 be arable land, waters, construction land, otherwise be the forest that vegetative coverage is high; When the NDVI value of the 9th phase is less than 0.25, be then waters, otherwise for ploughing and construction land; When the NDVI of the 2nd phase is less than 0.45, and when the EVI of the 9th phase is less than 0.32, be then construction land, otherwise for ploughing; 16th when the phase, NDVI was less than 0.45, be then spinney, otherwise be coniferous forest, broad-leaf forest, bamboo grove; When the NDVI of the 18th phase is less than 0.75, and when the EVI of the 16th phase is greater than 0.35, be then coniferous forest, otherwise be broad-leaf forest and bamboo grove; When the EVI of o. 11th is less than 0.65, and when the NDVI of the 15th phase is less than 0.75, be then bamboo grove, otherwise be then broad-leaf forest; Finally select the 2nd of NDVI the, 8,9,15,16, the 9th, 11,16 issues of 18 phases and EVI are according to synthesizing.

Described method, the step that described step 4 comprises is:

Opposite end tuple is divided and is carried out preliminary purification, minimal noise is carried out to the data obtained in step 3 and is separated conversion, pure index analysis, utilize N to tie up scatter diagram and carry out N dimension divergence analysis, tentatively determine terminal pixel, respective initial terminal pixel is obtained by tieing up in scatter diagram at N, the image after decision tree classification is updated to as region of interest, utilize decision tree classification to set up end member purification model to purify further to initial terminal pixel, the loose point of the inconsistent terminal pixel of the distribution plan obtained with decision tree classification is rejected or revised, obtains final terminal pixel.

Described method, the step that described step 5 comprises is:

Utilize final terminal pixel to ask for the spectral reflectivity of each atural object at each wave band, and in conjunction with GPS positioning sample place, each band spectrum reflectivity obtained by remote sensing image, determine the final atural object reflectance value of each atural object at each wave band.

Described method, the step that described step 6 comprises is:

Utilize each atural object at the final atural object reflectance value of each wave band, by linear spectral unmixing model

$R=\underset{k=1}{\overset{m}{\Σ}}{F}_k{r}_k+{\mathrm{\ϵ}}_k,$

Wherein m is end member number, r _kthe reflectivity of a kth end member group, F _kthe area ratio that in pixel, kth end-member composition is shared in pixel, ε _kfor the error of kth wave band, wherein the decomposition of the linear hybrid pixel of belt restraining will meet the ratio of component F of each pixel _ksum is 1, carries out Decomposition of Mixed Pixels to remote sensing image, obtains the abundance figure of each atural object.

Described method, the step that described step 7 comprises is:

Terminal unit Abundances is selected to be greater than the pixel of 0.5 often organizing in constitutional diagram picture, again often organizing constitutional diagram picture composition piece image, utilize maximum likelihood method to classify to classify to this image, obtain each Forest Types distribution plan, realize the information extraction of Decomposition of Mixed Pixels forest.

The technical solution adopted for the present invention to solve the technical problems is: carrying out on the pre-service bases such as projection transform, image mosaic, cutting, black removal band to the MODIS standardized product image adopted, first, utilize Forest Types phenology difference characteristic, adopt MOD13Q1 vegetation index product, by to normalized differential vegetation index NDVI and enhancement mode meta file EVI process, build NDVI, EVI time series section curve, by decision tree classification computer graphics sorting technique, build end member purification model; Secondly, adopt the MOD13Q1 vegetation index product building end member purification model to combine with MOD09A1 Reflectivity for Growing Season product and form new data, and the end member purification model of utilization structure carries out end member purifies further; Again, linear spectral unmixing model is adopted to decompose the new data be combined into; Finally, adopt MOD09A1 Reflectivity for Growing Season product and MOD13Q1 vegetation index product standard audio and video products, classified by maximum likelihood method, realize forest information extraction.The invention has the beneficial effects as follows, while the information of rapid extraction forest cover on a large scale, forest mapping precision and Forest Types accuracy of identification can be improved.MOD09A1 Reflectivity for Growing Season product and MOD13Q1 vegetation index product is only adopted in information extraction, satisfactory for result, simple.

Accompanying drawing explanation

Below in conjunction with drawings and Examples, the present invention is further described.

Fig. 1 is technical schematic diagram of the present invention;

Fig. 2 is the forest information extraction result figure of the embodiment of the present invention;

Fig. 3 is end member purification Decision-Tree Classifier Model figure of the present invention.

Embodiment

The present invention includes following steps:

Step one: data acquisition step, obtains forest remote sensing image;

Step 2: pre-treatment step, carries out pre-service to the remote sensing image that step one obtains, obtains pretreated remote sensing image data, pre-service comprises: carry out projection transform to remotely-sensed data, atmospheric correction, image mosaic, eliminate black stripe, geometry correction, study area is extracted;

Step 3: Forest Types phenology difference characteristic analytical procedure, finds out the vegetation index period obviously can distinguishing each Forest Types; Specifically, first GPS is utilized to obtain the coordinate information of (coniferous forest, broad-leaf forest, bamboo grove, spinney, waters, arable land, construction land), GIS software is utilized to obtain coniferous forest, broad-leaf forest, bamboo grove, spinney, waters, arable land, the gray-scale value of construction land field typical case's sampling point on remote sensing image, build vegetation index time series section curve, the vegetation index remotely-sensed data obviously can distinguishing each Forest Types found out by phenology analysis is combined with MOD09A1 Reflectivity for Growing Season product and forms new data.When the threshold value of the 8th phase NDVI be less than 0.6 be arable land, waters, construction land, otherwise be the forest that vegetative coverage is high; When the NDVI value of the 9th phase is less than 0.25, be then waters, otherwise for ploughing and construction land; When the NDVI of the 2nd phase is less than 0.45, and when the EVI of the 9th phase is less than 0.32, be then construction land, otherwise for ploughing; 16th when the phase, NDVI was less than 0.45, be then spinney, otherwise be coniferous forest, broad-leaf forest, bamboo grove; When the NDVI of the 18th phase is less than 0.75, and when the EVI of the 16th phase is greater than 0.35, be then coniferous forest, otherwise be broad-leaf forest and bamboo grove; When the EVI of o. 11th is less than 0.65, and when the NDVI of the 15th phase is less than 0.75, be then bamboo grove, otherwise be then broad-leaf forest; Therefore have selected the 2nd of NDVI the, 8,9,15,16, the 9th, 11,16 issues of 18 phases and EVI are according to synthesizing.

Step 4: end-member composition purification step, obtains terminal pixel component;

First end member purification model is built, utilize decision tree classification computer graphics sorting technique, according to the feature of each Forest Types phenology difference, utilize the vegetation index remotely-sensed data in period obviously can distinguishing each Forest Types found out, and utilize the vegetation index gray-scale value can distinguishing each atural object as the threshold value of Decision-Tree Classifier Model, obtain end member purification Decision-Tree Classifier Model as Fig. 3, in model, " ZYL " represents coniferous forest, " KYL " represents broad-leaf forest, " ZL " represents bamboo grove, " GML " represents spinney, " SY " represents waters, " GD " represents arable land, " JSYD " represents construction land, B2, B8, B9, B15, B16, B18 represents the 2nd of NDVI the respectively, 8, 9, 15, 16, 18 phases, b9, b11, b16, represent the 9th of EVI the respectively, 11, 16 phases.

Then opposite end tuple is divided and is carried out preliminary purification, minimal noise is carried out to the data obtained in step 3 and is separated conversion, pure index analysis, utilize N to tie up scatter diagram and carry out N dimension divergence analysis, tentatively determine terminal pixel, respective initial terminal pixel is obtained by tieing up in scatter diagram at N, the image after decision tree classification is updated to as region of interest, utilize decision tree classification to set up end member purification model to purify further to initial terminal pixel, the loose point that some do not meet the terminal pixel of rule and mistake (distribution plan namely obtained with decision tree classification is inconsistent) is rejected or revised, obtain final terminal pixel.

Step 5: end-member composition spectral reflectivity asks for step, obtains each reflectance; Namely utilize final terminal pixel to ask for the spectral reflectivity of each atural object at each wave band, and in conjunction with GPS positioning sample place, each band spectrum reflectivity obtained by remote sensing image, determine the final atural object reflectance value of each atural object at each wave band.

Step 6: linear spectral unmixing step, asks for each atural object distribution plan; Namely utilize each atural object at the final atural object reflectance value of each wave band, by linear spectral unmixing model wherein m is end member number, r _kthe reflectivity of a kth end member group, F _kthe area ratio that in pixel, kth end-member composition is shared in pixel, ε _kfor the error of kth wave band.Wherein the decomposition of the linear hybrid pixel of belt restraining will meet the ratio of component F of each pixel _ksum is 1, carries out Decomposition of Mixed Pixels to remote sensing image, obtains the abundance figure of each atural object.

Step 7: forest information extracting step, obtains each atural object distribution plan.Namely terminal unit Abundances is selected to be greater than the pixel of 0.5 often organizing in constitutional diagram picture, again often organizing constitutional diagram picture composition piece image, utilize maximum likelihood method to classify to classify to this image, obtain each Forest Types distribution plan, realize the information extraction of Decomposition of Mixed Pixels forest.

See Fig. 1, when implementing, need to carry out the pre-service such as projection transform, image mosaic, cutting, black removal band to MOD09A1 Reflectivity for Growing Season product and MOD13Q1 vegetation index product, by to MOD13Q1 normalized differential vegetation index NDVI and strengthen vegetation index EVI, build NDVI, EVI time series section curve, analyze the threshold value drawing and build decision tree purification model.And to combine with MOD09A1 Reflectivity for Growing Season product produces new data building the threshold data product of decision-tree model, utilize linear spectral unmixing model to decompose, recycle maximum likelihood method and carry out classification extraction forest information.

In the embodiment depicted in figure 2, choosing Hunan Province is object, adopts the MOD09A1 Reflectivity for Growing Season product of 2009 and MOD13Q1 vegetation index product to carry out the pre-service such as projection transform, image mosaic, cutting, black removal band.First, utilize Forest Types phenology difference characteristic, adopt MOD13Q1 vegetation index product, by processing normalized differential vegetation index NDVI and enhancement mode meta file EVI (annual each 23 phases), build NDVI, EVI time series section curve, filter out the data of optimum decision tree disaggregated model, build end member purification model; Secondly, adopt the MOD13Q1 vegetation index product building end member purification model to combine with MOD09A1 Reflectivity for Growing Season product and form new data, and carry out minimum noise separation (MNF), pixel purity index calculates (PPI), density slice, and N ties up scatter diagram, end-member composition preliminary purification process, and utilize decision tree purification model to purify further, build the spectral reflectivity marking each wave band in conjunction with image, draw final end-member composition reflectivity; Again, linear spectral unmixing model is adopted to decompose the data be combined into; Finally, the data adopting MOD09A1 Reflectivity for Growing Season product to become with the MOD13Q1 vegetation index product mix screened, are classified by maximum likelihood method, realize forest information extraction, and process exports Hunan Province's forest distribution map.

The MODIS Decomposition of Mixed Pixels forest information extracting method that the present invention proposes has increased substantially forest mapping and Forest Types accuracy of identification, and classification results is desirable, has important practical dissemination.

Claims (9)

1. a MODIS Decomposition of Mixed Pixels forest information extracting method, is characterized in that, comprise the steps:

Step one: data acquisition step, obtains forest remote sensing image;

Step 2: pre-treatment step, carries out pre-service to the remote sensing image that step one obtains, obtains pretreated remote sensing image data;

Step 3: Forest Types phenology difference characteristic analytical procedure, finds out the vegetation index period obviously can distinguishing each Forest Types;

Step 4: end-member composition purification step, obtains terminal pixel component;

Step 5: end-member composition spectral reflectivity asks for step, obtains each reflectance;

Step 6: linear spectral unmixing step, asks for each atural object distribution plan;

Step 7: forest information extracting step, obtains each atural object distribution plan.

2. method according to claim 1, is characterized in that, in described step 2, the pre-service of remotely-sensed data comprises: carry out projection transform to remotely-sensed data, atmospheric correction, image mosaic, eliminates black stripe, geometry correction, and study area is extracted.

3. method according to claim 2, is characterized in that, the step that described step 3 comprises is:

GPS is utilized to obtain the coordinate information of typical sampling point, GIS software is utilized to obtain the gray-scale value of typical sampling point on remote sensing image, build vegetation index time series section curve, the data that the vegetation index remotely-sensed data obviously can distinguishing each Forest Types then found out by phenology analysis and MOD09A1 Reflectivity for Growing Season product combined together are obtained.

4. method according to claim 3, is characterized in that, described typical sampling point comprises coniferous forest, broad-leaf forest, bamboo grove, spinney, waters, arable land and construction land.

5. method according to claim 3, is characterized in that, the data that the described vegetation index remotely-sensed data obviously can distinguishing each Forest Types found out by phenology analysis and MOD09A1 Reflectivity for Growing Season product combined together are obtained comprise the following steps:

First end member purification model is built, utilize decision tree classification computer graphics sorting technique, according to the feature of each Forest Types phenology difference, utilize the vegetation index remotely-sensed data in period obviously can distinguishing each Forest Types found out, and utilize the vegetation index gray-scale value can distinguishing each atural object as the threshold value of Decision-Tree Classifier Model: namely

When the threshold value of the 8th phase NDVI be less than 0.6 be arable land, waters, construction land, otherwise be the forest that vegetative coverage is high; When the NDVI value of the 9th phase is less than 0.25, be then waters, otherwise for ploughing and construction land; When the NDVI of the 2nd phase is less than 0.45, and when the EVI of the 9th phase is less than 0.32, be then construction land, otherwise for ploughing; 16th when the phase, NDVI was less than 0.45, be then spinney, otherwise be coniferous forest, broad-leaf forest, bamboo grove; When the NDVI of the 18th phase is less than 0.75, and when the EVI of the 16th phase is greater than 0.35, be then coniferous forest, otherwise be broad-leaf forest and bamboo grove; When the EVI of o. 11th is less than 0.65, and when the NDVI of the 15th phase is less than 0.75, be then bamboo grove, otherwise be then broad-leaf forest; Finally select the 2nd of NDVI the, 8,9,15,16, the 9th, 11,16 issues of 18 phases and EVI are according to synthesizing.

6. method according to claim 3, is characterized in that, the step that described step 4 comprises is:

Opposite end tuple is divided and is carried out preliminary purification, minimal noise is carried out to the data obtained in step 3 and is separated conversion, pure index analysis, utilize N to tie up scatter diagram and carry out N dimension divergence analysis, tentatively determine terminal pixel, respective initial terminal pixel is obtained by tieing up in scatter diagram at N, the image after decision tree classification is updated to as region of interest, utilize decision tree classification to set up end member purification model to purify further to initial terminal pixel, the loose point of the inconsistent terminal pixel of the distribution plan obtained with decision tree classification is rejected or revised, obtains final terminal pixel.

7. method according to claim 6, is characterized in that, the step that described step 5 comprises is:

Utilize final terminal pixel to ask for the spectral reflectivity of each atural object at each wave band, and in conjunction with GPS positioning sample place, each band spectrum reflectivity obtained by remote sensing image, determine the final atural object reflectance value of each atural object at each wave band.

8. method according to claim 7, is characterized in that, the step that described step 6 comprises is:

Utilize each atural object at the final atural object reflectance value of each wave band, by linear spectral unmixing model

$R=\underset{k=1}{\overset{m}{\mathrm{\Σ}}}{F}_k{r}_k+{\mathrm{\ϵ}}_k,$

Wherein m is end member number, r _kthe reflectivity of a kth end member group, F _kthe area ratio that in pixel, kth end-member composition is shared in pixel, ε _kfor the error of kth wave band, wherein the decomposition of the linear hybrid pixel of belt restraining will meet the ratio of component F of each pixel _ksum is 1, carries out Decomposition of Mixed Pixels to remote sensing image, obtains the abundance figure of each atural object.

9. method according to claim 8, is characterized in that, the step that described step 7 comprises is:

Terminal unit Abundances is selected to be greater than the pixel of 0.5 often organizing in constitutional diagram picture, again often organizing constitutional diagram picture composition piece image, utilize maximum likelihood method to classify to classify to this image, obtain each Forest Types distribution plan, realize the information extraction of Decomposition of Mixed Pixels forest.