CN113538388B - Arable land loss assessment method based on MODIS NDVI time sequence data - Google Patents
Arable land loss assessment method based on MODIS NDVI time sequence data Download PDFInfo
- Publication number
- CN113538388B CN113538388B CN202110837952.9A CN202110837952A CN113538388B CN 113538388 B CN113538388 B CN 113538388B CN 202110837952 A CN202110837952 A CN 202110837952A CN 113538388 B CN113538388 B CN 113538388B
- Authority
- CN
- China
- Prior art keywords
- land
- time sequence
- breakpoint
- sequence data
- formula
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Active
Links
Images
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T7/00—Image analysis
- G06T7/0002—Inspection of images, e.g. flaw detection
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T7/00—Image analysis
- G06T7/60—Analysis of geometric attributes
- G06T7/62—Analysis of geometric attributes of area, perimeter, diameter or volume
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T2207/00—Indexing scheme for image analysis or image enhancement
- G06T2207/10—Image acquisition modality
- G06T2207/10032—Satellite or aerial image; Remote sensing
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T2207/00—Indexing scheme for image analysis or image enhancement
- G06T2207/30—Subject of image; Context of image processing
- G06T2207/30181—Earth observation
- G06T2207/30184—Infrastructure
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T2207/00—Indexing scheme for image analysis or image enhancement
- G06T2207/30—Subject of image; Context of image processing
- G06T2207/30181—Earth observation
- G06T2207/30188—Vegetation; Agriculture
Abstract
The invention discloses a farmland loss evaluation method based on MODISNDVI time sequence data, and belongs to the technical field of remote sensing image processing. Firstly, extracting pixels which are converted from cultivated land into construction land by using the existing ground surface covering product; and then, an NDVI time sequence is constructed by utilizing MODISNDVI data synthesized in the month, and the NDVI data is decomposed by utilizing a mathematical model to obtain trend components, seasonal components and break points existing in the middle, wherein the positions of the break points represent the time when the land type changes. The method can obtain the time for converting the farmland into the construction land, and accordingly, the area of the corresponding region from the farmland into the construction land is estimated every year.
Description
Technical Field
The invention belongs to the field of remote sensing image processing and remote sensing change detection, and particularly relates to a cultivated land loss evaluation method based on MODIS NDVI time sequence data, which is used for detecting the cultivated land loss change caused by city expansion by utilizing MODIS NDVI time sequence data and acquiring the specific time for converting the cultivated land into a construction land.
Background
In 2001-2005, the agricultural land occupied by urban expansion reaches 2178km2 per year, and most of the occupied land is high-quality cultivated land around plain areas or residential areas; the Chinese development report 2010 predicts that the Chinese population scale reaches the peak and the urbanization rate reaches 65% in 20-30 years in the future. The rapid urbanization process further reduces the agricultural land use to a large extent, and forms a potential threat to the food safety of China. The method can accurately and objectively acquire the agricultural land loss condition caused by urban expansion every year, and can provide effective scientific basis for making policies such as urban development and land management. However, the conventional method for counting the current situation of agricultural land is land survey, which consumes a lot of financial resources and manpower and cannot be frequently carried out; meanwhile, the land survey duration period is long, and the data have different degrees of hysteresis. Remote sensing has advantages such as large-scale, short cycle, continuity, can make up the shortcoming of land investigation well.
Currently, the method for analyzing the farmland change based on the remote sensing image is mainly a change detection method, and many researchers provide a remote sensing image change detection model, which is summarized as follows:
1) The basic land utilization change detection method is a multi-temporal remote sensing image classification method, namely, multi-period remote sensing images are subjected to earth surface coverage classification, and then change areas are extracted, wherein the classification method comprises a support vector machine, a maximum likelihood classification method, a deep learning method and the like, however, the accuracy of the current classification algorithm is still not high due to the complex scene of the remote sensing images, and a single classification model cannot be suitable for all areas;
2) The second change detection method is based on Principal Component Analysis (PCA) of multi-temporal remote sensing images to reduce the influence caused by classification errors, the method does not depend on image classification, but is influenced by the climate, the features of the same earth surface covered on the images change along with the season, and therefore the method is time-efficient;
3) Currently, with the development of deep learning, many researchers have proposed a remote sensing image change detection method based on deep learning, such as a change detection method based on ResNet, which is limited by the number of samples, cannot be applied to all regions, and has limited accuracy.
The change detection method can only know whether the land type of the farmland changes, and cannot determine the specific time when the farmland changes.
Disclosure of Invention
The invention aims to solve the technical problem of extracting the spatial and temporal change condition of farmland loss caused by urban expansion, and provides an farmland loss evaluation method based on MODIS NDVI time sequence data.
The technical scheme adopted by the invention is as follows:
a farmland loss assessment method based on MODIS NDVI time sequence data comprises the following steps:
(1) Obtaining MODIS NDVI time sequence data, respectively converting and projecting the time sequence data to a geographic coordinate system, and superposing the time sequence data together according to dates to generate time sequence data;
(2) Extracting pixels of which the first-stage land covering type is cultivated land and the later-stage land covering type is construction land by using the existing land utilization product;
(3) Fitting the pixels extracted in the step (2) pixel by using a time sequence model;
(4) Judging whether a breakpoint exists in MODIS NDVI time sequence data and the breakpoint time by using an AIC (automatic interface component) criterion according to a fitting result; and if the breakpoint exists, judging whether the earth surface coverage type is changed from cultivated land to construction land.
Further, the time sequence model in step (3) is:
Y v ＝T v +S v +R v (1)
wherein:
T v ＝α+βt (3)
in the formula, T v A representative trend component, expressed by a linear function; s v Representing a seasonal component, expressed by a Fourier series of trigonometric functions; r v Representing residual error, alpha, beta being coefficients of linear functions, gamma k And theta k F is the number of MODIS NDVI observations per year, k is the coefficient determining the period of the trigonometric function, and t is an independent variable representing the t-th observation.
Further, the fitting mode in the step (3) is specifically as follows:
if the earth surface coverage does not change, the MODIS NDVI time sequence data has no breakpoint, so that:
obtaining residual square sum RSS1 according to the fitting result;
if groundWhen the table coverage changes, 1 breakpoint exists in the MODIS NDVI time sequence data, and the breakpoint position is set as t * ，f＜t * < L-f, so:
in the formula, f is a period, L is the total amount of MODIS NDVI time sequence data observation, alpha 1 、β 1 、α 2 、β 2 Respectively coefficients of two linear functions, gamma i 、θ i 、γ j 、θ j For coefficients, i, j are coefficients determining the period of the trigonometric function before and after the break point, t * The position of the breakpoint is shown, t is an independent variable and represents the t-th observed value;
traversal of t using equation (5) * The minimum value RSS2 of the sum of the squares of the residuals is determined, and t at this time is recorded * Value of (a)
I.e. the breakpoint time.
Further, in the step (4), the AIC criterion is used to determine whether a breakpoint exists in the MODIS NDVI timing data and a breakpoint time, specifically:
the AIC criterion expression is as follows:
AIC＝2l+Lln(RSS/L) (6)
in the formula, l is the number of parameters in the formula (4) or the formula (5), the number of parameters in the formula (4) or the formula (5) is substituted into the formula (6) to respectively obtain two AIC values AIC1 and AIC2, and if the AIC1 is adopted, the two AIC values AIC1 and AIC2 are respectively obtained>AIC2, a breakpoint exists in the time sequence data, and the breakpoint time is
The breakpoint moment is the moment when the earth surface coverage type changes.
Further, the step (4) of judging whether the surface covering type is changed from the cultivated land to the construction land or not specifically comprises the following steps:
if the ground cover type is changed from arable land to construction land, the amplitude of the seasonal component before and after the breakpoint is reduced, i.e.:
if equation (7) is satisfied, the ground surface coverage type is changed from the cultivated land to the construction land.
Compared with the background technology, the invention has the following beneficial effects:
the method can obtain the time for each pixel to be converted from the cultivated land to the construction land, and accordingly, the area for each pixel to be converted from the cultivated land to the construction land is counted.
Drawings
FIG. 1 is a flow chart of the tillage loss evaluation according to the embodiment of the invention.
Fig. 2 is a schematic diagram illustrating decomposition of MODIS NDVI timing data according to an embodiment of the present invention.
FIG. 3 shows the area of year-by-year cultivation lost in Shandong province in accordance with an embodiment of the present invention.
Detailed Description
The following description of the embodiments of the present invention is provided in order to better understand the present invention for those skilled in the art with reference to the accompanying drawings. It is to be expressly noted that in the following description, a detailed description of known functions and designs will be omitted when it may obscure the subject matter of the present invention.
Specifically, the method comprises the following steps:
(1) And constructing the MODIS NDVI time sequence data of the research area. Downloading MODIS NDVI data with the time resolution higher than 1 month, respectively converting and projecting the MODIS NDVI data to a geographic coordinate system, and then superposing all the MODIS NDVI data together according to dates to generate time sequence data; and cutting out MODIS NDVI time sequence data of the research area by using the vector data of the research area.
(2) Extracting pixels of which the first-stage land covering type is cultivated land and the later-stage land covering type is construction land by utilizing the existing land utilization product;
(3) For the pixels extracted in (2), modeling and fitting are performed pixel by pixel according to the following formula:
Y v ＝T v +S v +R v (1)
wherein, T v A representative trend component, expressed by a linear function; s v Representing a seasonal component, expressed by a Fourier series of trigonometric functions; r is v Representing the residual.
Wherein:
T v ＝α+βt (3)
where α, β are coefficients of linear functions, γ k And theta k F is the number of MODIS NDVI observations per year, k is the coefficient determining the period of the trigonometric function, and t is an independent variable representing the t-th observation.
1) If the earth surface coverage does not change, the MODIS NDVI time sequence data has no breakpoint, so that:
the fitting results in a residual sum of squares RSS1.
2) If the earth surface coverage changes, the earth surface coverage changes are not frequent according to general conditions, 1 breakpoint exists in the MODIS NDVI time sequence data, and the breakpoint position is assumed to be t * (f＜t * < L-f) (generally, f takes the value of one period; l is the total number of observations of MODIS NDVI timing data), so:
in the formula, f is a period, L is the total number of MODIS NDVI time sequence data observation, alpha 1 、β 1 、α 2 、β 2 Respectively coefficients of two linear functions, gamma i 、θ i 、γ j 、θ j Is a coefficient, i, j are before and after the breakpointDetermining the coefficient of the period of the trigonometric function, t * The position of the breakpoint is shown, t is an independent variable and represents the t-th observed value;
traversal of t using equation (5) * To find the residual sum of squares minimum RSS2 and record t at this time * Is taken as
(4) And judging whether a breakpoint exists in the MODISNDVI time sequence by using the AIC criterion.
The AIC criteria expression is as follows:
AIC＝2l+Lln(RSS/L) (6)
wherein l is the number of parameters of formula (4) (no breakpoint exists in MODIS NDVI timing data) or formula (5) (one breakpoint exists in MODIS NDVI timing).
Substituting the parameter numbers of RSS1 and formula (4) into formula (6) to obtain AIC1, and substituting the parameter numbers of RSS2 and formula (5) into formula (6) to obtain AIC2; if AIC1<AIC2, there is no breakpoint in the time sequence data, otherwise, if AIC1>AIC2, a breakpoint exists in the time sequence data, and the breakpoint time is
The time is the time when the earth's surface coverage type changes.
Determining whether the surface coverage change belongs to a change from arable land to a building. If the ground cover type is changed from arable land to construction land, the amplitude of the seasonal component before and after the breakpoint is reduced, i.e.:
if the surface coverage change is to be transferred from arable land to a structure, the mask map is set to 1, otherwise it is set to 0.
The following is a more specific example:
as shown in FIG. 1, a farmland loss evaluation model based on MODIS time series data comprises the following steps:
1. constructing MODISNDVI time sequence data in research area
The method comprises the steps of downloading monthly synthetic MODISNDVI data with row and column numbers of h27v05 from 1 month to 2009 in 2001, wherein the annual observation quantity is 23, 207 observation values are totally obtained, and the spatial resolution is 250 meters; and converting the projection into a geographic coordinate system, then superposing the projection according to time, and clipping according to the administrative division of Shandong province.
2. And extracting pixels with the land cover type of arable land in 2000 and the land cover type of construction land in 2010 by using products of GlobeLand30-2000 and GlobeLand30-2010, constructing a land use change mask map, and resampling the mask map to 250 meters by using a mode polymerization method.
3. Pixel by pixel modeling and fitting was done according to equation (1) where f =23. Fig. 2 gives the fitting result for a certain pixel.
4. And (3) judging whether a breakpoint exists in the time sequence data or not by utilizing an AIC (automatic information center) criterion in the modeling and fitting process of each pixel, and determining a final fitting mode. In fig. 2, the breakpoint is determined to exist according to the AIC criterion, and the breakpoint position is the 94 th observation value, which is month 2 of 2005.
5. And judging the seasonal component amplitude change before and after the breakpoint for each pixel with changed earth surface coverage, and if the amplitude is reduced, confirming that the earth surface coverage of the pixel is changed from cultivated land to construction land. The amplitude of the seasonal component in fig. 2 is reduced after the surface coverage is changed, and it can be determined that the arable land is changed into the construction land; namely, the pixel is changed from farmland to construction land in 2 months of 2005.
6. The area of the land for construction changed from the cultivated land in Shandong province is counted according to the year, as shown in figure 3.
In a word, the invention utilizes MODISNDVI time sequence data to extract the pixels converted from the cultivated land into the construction land and obtain the time for converting the cultivated land into the construction land.
Although the illustrative embodiments of the present invention have been described in order to facilitate those skilled in the art to understand the present invention, it is to be understood that the present invention is not limited to the scope of the embodiments, and that various changes may be made apparent to those skilled in the art as long as they are within the spirit and scope of the present invention as defined and defined in the appended claims, and all matters of the invention using the inventive concepts are protected.
Claims (2)
1. A farmland loss assessment method based on MODIS NDVI time sequence data is characterized by comprising the following steps:
(1) Obtaining MODIS NDVI time sequence data, respectively converting and projecting the time sequence data to a geographic coordinate system, and superposing the time sequence data together according to dates to generate time sequence data;
(2) Extracting pixels of which the first-stage land covering type is cultivated land and the later-stage land covering type is construction land by using the existing land utilization product;
(3) Fitting the pixels extracted in the step (2) pixel by using a time sequence model; wherein, the time sequence model is as follows:
Y v ＝T v +S v +R v (1)
T v ＝α+βt(3)
in the formula, T v A representative trend component, expressed by a linear function; s v Representing a seasonal component, expressed by a Fourier series of trigonometric functions; r v Representing residual error, alpha, beta being coefficients of linear functions, gamma k And theta k F is the number of MODIS NDVI observations per year, k is the coefficient for determining the period of the trigonometric function, and t is an independent variable and represents the t-th observation value;
the fitting mode is specifically as follows:
if the earth surface coverage does not change, no breakpoint exists in MODIS NDVI time sequence data, so that:
obtaining residual square sum RSS1 from the fitting result;
if the earth surface coverage changes, 1 breakpoint exists in MODIS NDVI time sequence data, and the breakpoint position is set as t * ，f＜t * < L-f, so:
in the formula, f is a period, L is the total number of MODIS NDVI time sequence data observation, alpha 1 、β 1 、α 2 、β 2 Respectively coefficients of two linear functions, gamma i 、θ i 、γ j 、θ j For coefficients, i, j are coefficients determining the period of the trigonometric function before and after the break point, t * The position of the breakpoint is shown, t is an independent variable and represents the t-th observed value;
traversal of t using equation (5) * The minimum value RSS2 of the sum of the squares of the residuals is determined, and t at this time is recorded * Value of
Namely the breakpoint moment;
(4) Judging whether a breakpoint exists in MODIS NDVI time sequence data and the breakpoint time by using an AIC (automatic interface component) criterion according to a fitting result; if the breakpoint exists, judging whether the earth surface coverage type is changed from the cultivated land to the land for building; the method for judging whether the MODIS NDVI time sequence data has the breakpoint and the breakpoint time by using the AIC criterion specifically comprises the following steps of:
the AIC criteria expression is as follows:
AIC＝2l+Lln(RSS/L)(6)
in the formula, l is the number of parameters in the formula (4) or the formula (5), the number of parameters in the formula (4) or the formula (5) is substituted into the formula (6) to respectively obtain two AIC values AIC1 and AIC2, and if the AIC1 is adopted, the two AIC values AIC1 and AIC2 are respectively obtained>AIC2, a breakpoint exists in the time sequence data, and the breakpoint time is
At the moment of the breakpointThe time when the type of surface coverage changes.
2. The method for assessing the arable land loss based on MODIS NDVI time series data according to claim 1, wherein the step (4) is performed to determine whether the surface coverage type is changed from arable land to a construction land, specifically:
if the ground cover type is changed from arable land to construction land, the amplitude of the seasonal component before and after the breakpoint is reduced, i.e.:
if equation (7) is satisfied, the ground surface coverage type is changed from the cultivated land to the construction land.
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202110837952.9A CN113538388B (en) | 2021-07-23 | 2021-07-23 | Arable land loss assessment method based on MODIS NDVI time sequence data |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202110837952.9A CN113538388B (en) | 2021-07-23 | 2021-07-23 | Arable land loss assessment method based on MODIS NDVI time sequence data |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN113538388A CN113538388A (en) | 2021-10-22 |
| CN113538388B true CN113538388B (en) | 2022-10-11 |
Family
ID=78089411
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202110837952.9A Active CN113538388B (en) | 2021-07-23 | 2021-07-23 | Arable land loss assessment method based on MODIS NDVI time sequence data |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN113538388B (en) |
Citations (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN111768101A (en) * | 2020-06-29 | 2020-10-13 | 北京师范大学 | Remote sensing farmland change detection method and system considering phenological characteristics |
Family Cites Families (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN104318270A (en) * | 2014-11-21 | 2015-01-28 | 东北林业大学 | Land cover classification method based on MODIS time series data |
| CN113607284B (en) * | 2017-01-23 | 2022-06-24 | 北京师范大学 | Method for distinguishing potential fire point by using BFAST algorithm |
| CN108564002B (en) * | 2018-03-22 | 2020-08-21 | 中国科学院遥感与数字地球研究所 | Method and system for detecting time sequence change of remote sensing image |
| CN112001809A (en) * | 2020-07-31 | 2020-11-27 | 中科海慧（天津）科技有限公司 | Method for acquiring farmland returning information of agriculture and forestry area |
-
2021
- 2021-07-23 CN CN202110837952.9A patent/CN113538388B/en active Active
Patent Citations (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN111768101A (en) * | 2020-06-29 | 2020-10-13 | 北京师范大学 | Remote sensing farmland change detection method and system considering phenological characteristics |
Also Published As
| Publication number | Publication date |
|---|---|
| CN113538388A (en) | 2021-10-22 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| Xue et al. | Phenology-driven land cover classification and trend analysis based on long-term remote sensing image series | |
| Zhang et al. | Coastal wetland vegetation classification with a Landsat Thematic Mapper image | |
| CN103914678B (en) | Abandoned land remote sensing recognition method based on texture and vegetation indexes | |
| Camus et al. | Regional analysis of multivariate compound coastal flooding potential around Europe and environs: sensitivity analysis and spatial patterns | |
| Soria-Ruiz et al. | Land-cover classification using radar and optical images: a case study in Central Mexico | |
| CN111982822B (en) | Long-time sequence high-precision vegetation index improvement algorithm | |
| Zhang et al. | Spatial domain bridge transfer: An automated paddy rice mapping method with no training data required and decreased image inputs for the large cloudy area | |
| van Oort et al. | Spatial variability in classification accuracy of agricultural crops in the Dutch national land-cover database | |
| Yang et al. | Spatial-temporal dynamic monitoring of vegetation recovery after the Wenchuan earthquake | |
| CN112329265A (en) | Satellite remote sensing rainfall refinement space estimation method and system | |
| Patidar et al. | A rule-based spectral unmixing algorithm for extracting annual time series of sub-pixel impervious surface fraction | |
| CN114120137B (en) | Time-sequence vegetation remote sensing image-based wetland element time-space evolution monitoring method | |
| Zhang | A time-series approach to detect urbanized areas using biophysical indicators and landsat satellite imagery | |
| CN113901348A (en) | Oncomelania snail distribution influence factor identification and prediction method based on mathematical model | |
| CN113538388B (en) | Arable land loss assessment method based on MODIS NDVI time sequence data | |
| Faid et al. | Monitoring land-use change-associated land development using multitemporal Landsat data and geoinformatics in Kom Ombo area, South Egypt | |
| Zhou et al. | Modelling spatio-temporal pattern of landuse change using multi-temporal remotely sensed imagery | |
| CN115797501A (en) | Time-series forest age mapping method combining forest disturbance and recovery events | |
| Olmanson et al. | Use of Landsat imagery to develop a water quality atlas of Minnesota’s 10,000 lakes | |
| Cheng et al. | Technical framework of feature extraction based on pixel-level SAR image time series | |
| CN112052720B (en) | High-space-time normalization vegetation index NDVI fusion model based on histogram clustering | |
| Zhu et al. | A coupled temporal-spectral-spatial multidimensional information change detection framework method: a case of the 1990-2020 Tianjin, China | |
| Yang et al. | Remote sensing and geospatial analysis for landscape pattern characterization | |
| Liu et al. | Fish-pond change detection based on short term time series of RADARSAT images and object-oriented method | |
| CN113609899B (en) | Remote sensing time sequence analysis-based tilling land information positioning display system |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant |