A large-scale dynamic impervious surface coverage estimation method based on time series optimization

Abstract

The present invention discloses a large-scale dynamic impervious surface coverage estimation method based on time series optimization, which comprises the following steps: S1: multi-source remote sensing data acquisition and preprocessing in the monitoring area; S2: time series impervious surface coverage sample preparation; S3: impervious surface coverage feature set construction; S4: impervious surface coverage estimation model construction and spatiotemporal information extraction; S5: impervious surface coverage time series logic optimization; S6: result verification. The method comprehensively utilizes multi-source remote sensing data to prepare time series dynamic impervious surface samples, obtains long-term series impervious surface coverage results through feature sets constructed by multi-source information, and corrects dynamic impervious surface changes based on time series logic optimization. It is a process for extracting impervious surface information with strong operability, high accuracy, stability and reliability, and has a positive contribution to large-scale impervious surface inversion.

Description

Technical Field

[0001] The invention relates to remote sensing information interpretation and remote sensing image processing technology, and specifically refers to a large-scale dynamic impervious surface coverage estimation method based on time series optimization. Background Art

[0002] Impervious surface is a kind of surface covering material that prevents water from penetrating into the soil, mainly including buildings, roads, parking lots, etc. Impervious surface extraction refers to the identification and extraction of artificial surfaces in cities or other areas from satellite images, aerial photos or ground data through remote sensing or other technical means. The increase of impervious surface will cause a series of environmental and urban problems, such as urban heat island, urban flooding, regional rainfall, carbon emissions, etc. Specifically, the expansion of impervious surface directly affects the urban heat island effect and changes the local climate environment. A higher proportion of impervious surface will lead to temperature rise, which will increase energy consumption and affect air quality. With the increase of impervious surface area, rainwater cannot penetrate effectively, which is very likely to form waterlogging and flood disasters. The expansion of impervious surface not only aggravates local environmental problems, but also may have a chain reaction on global climate change by affecting water vapor circulation and carbon emissions. At the same time, the dynamic change information of impervious surface is also an important reliance for describing the scale and development path of cities. The distribution and proportion of impervious surface are important indicators of urban land use and reflect the urbanization process. Therefore, monitoring impervious surfaces can help urban planners develop more sustainable urban development strategies and rationally allocate land resources. Through impervious surface extraction, it can help governments and relevant departments develop targeted environmental protection policies and urban development plans.

[0003] Therefore, the extraction of impervious surface information has far-reaching significance for urban management and ecological protection, and it is very important to draw and master accurate and reasonable time-series dynamic impervious surfaces.

[0004] Large-scale dynamic impervious surface changes are limited by the development of data and technology. Problems such as missing data and small coverage will affect the mapping of long-term impervious surfaces, and the low adaptability of methods and technologies will reduce the feasibility of large-scale impervious surface mapping.

[0005] At present, there are two main methods for extracting time-series impervious surfaces:

[0006] The first category is to obtain impervious surface information based on time series change detection methods, which usually refers to the use of long-term series remote sensing images such as Landsat, MODIS, and Sentinel to invert impervious surfaces based on surface cover changes reflected by spectral information. Common methods include mutation point methods for separating trend and seasonal terms, continuous change detection, and classification. These methods collect long-term continuously distributed remote sensing images, identify the time of possible changes through spectral information, and further extract time-series impervious surfaces. However, methods based on time series change detection require the acquisition of remote sensing data that is densely distributed in time, and the data quality will also affect the accuracy of the extraction results.

[0007] The second category is to obtain impervious surface information based on multi-temporal independent extraction methods, which usually refers to the use of long-term remote sensing images, high-resolution remote sensing images, etc. to invert the spatial distribution of impervious surfaces at each independent time. Commonly used methods include spectral mixture decomposition, index extraction, traditional image classification, machine learning, and deep learning. These methods extract impervious surfaces in each independent period. Machine learning and deep learning models are feasible methods, however, the quality of samples affects the prediction effect, and low-quality samples will lead to unreasonable results in the model. Early impervious surface prediction models relied on artificially produced samples as input samples, which required huge manpower and material resources. The method of deriving impervious surface samples using existing time series data can improve the accuracy of sample production, but there are differences in existing data, and the accuracy and reliability of these data need to be considered, and then appropriate evaluation and screening are carried out. In addition, for the time series impervious surface coverage, the required samples are continuous variables, which are affected by more time and space errors.

[0008] The method based on multi-temporal independent extraction needs to take into account the temporal changes of impervious surfaces and correct the temporal changes of each pixel. The temporal optimization of single impervious surface coverage refers to the phenomenon that the impervious surface remains unchanged or increases. This method is based on the assumption that the impervious surface will not change inversely. However, the changes in the impervious surface are complex.

[0009] In summary, impervious surface is an important indicator for detecting ecological environment and urban development, reflecting the regional ecological quality and urban pattern changes. Reasonable data and technology are more important for extracting impervious surface information.

[0010] It is still difficult to extract large-scale time series dynamic impervious surface information, and the existing methods have certain shortcomings: data with small coverage and short duration cannot realize large-scale and long-time series impervious surface monitoring; the original quality of the data will affect the accuracy of methods based on time series change detection; methods based on multi-phase independent extraction require time series correction; a single impervious surface coverage time series optimization method cannot express the characteristics of impervious surface changes under complex conditions.

[0011] Therefore, a large-scale dynamic impervious surface coverage estimation method based on time series optimization is proposed to solve the above problems. Summary of the invention

[0012] The purpose of this application is to provide a large-scale dynamic impervious surface coverage estimation method based on time series optimization to solve the problems raised in the above background technology. This method is particularly suitable for mapping the spatiotemporal distribution of impervious surfaces with long time series and large coverage.

[0013] The present application provides a large-scale dynamic impervious surface coverage estimation method based on time series optimization, which is characterized by comprising the following steps:

[0014] S1: Acquisition and preprocessing of multi-source remote sensing data in the monitoring area: obtain long-term Landsat data and MODIS data in the monitoring area and remove outliers, obtain time-series night light NTL data and multi-source impervious surface product data;

[0015] S2: Production of time series impervious surface coverage samples: Obtain time series impervious surface coverage samples through two methods: exporting impervious surface products and supplementing time series samples;

[0016] S3: Construction of impervious surface coverage feature set: Based on the time series MODIS data and time series nighttime light NTL data, the impervious surface coverage estimation feature set is constructed, including seasonal features, statistical features, auxiliary features and index features;

[0017] S4: Construction of impervious surface coverage estimation model and extraction of spatiotemporal information: Input the temporal impervious surface coverage samples and the corresponding impervious surface coverage feature set into the impervious surface coverage estimation model for training, and estimate the temporal impervious surface coverage results;

[0018] S5: Temporal logic optimization of impervious surface coverage: carry out temporal change correction according to the temporal impervious surface coverage under four impervious surface coverage change trends, the four impervious surface coverage change trends include: increase, decrease, unchanged and others;

[0019] S6: Result verification: Verify the accuracy of the extraction results based on the accuracy assessment method.

[0020] Preferably, in step S1, remote sensing images of the monitoring area in each year are collected through the GEE platform, including Landsat surface reflectance data Landsat5 / 7 / 8 and MODIS surface reflectance data MOD09A1, to form a long-term remote sensing image sequence, and outliers are detected pixel by pixel according to the quality assessment bands of Landsat and MODIS, and abnormal pixels are masked;

[0021] The multi-source impervious surface products are obtained through the GEE platform and the Zenodo platform, including: global artificial impervious surface annual dynamic data products, global impervious surface products, global impervious surface dynamic data sets, global annual urban dynamic data, and global human settlement layer products. The multi-source impervious surface products are converted into impervious surface coverage consistent with the time span and spatial resolution of remote sensing images. The calculation formula for impervious surface coverage C is as follows:

[0022]

[0023] Among them, sum(AS) is the area of ​​the impervious surface in the grid, and length is the length of the grid, that is, the length of the pixel side.

[0024] Preferably, in step S21, a temporally stable impervious surface coverage sample is obtained from multi-source impervious surface products by a threshold method. If the temporal fluctuation range of the impervious surface coverage of the same product at the pixel scale and the difference in impervious surface coverage between different products are both within a preset threshold range, the pixel is identified as a stable pixel as a temporally stable impervious surface coverage sample.

[0025] Preferably, in step S3, the seasonal features include summer features and winter features, which are synthesized by the time series MOD09A1 through the time series fusion method, specifically including the red band B1, the near infrared band B2, the blue band B3, the green band B4, the near infrared band B5, the near infrared band B6, and the short-wave infrared band B7. The MOD09A1 images of summer and winter are respectively selected to generate the spectral time series synthesis features of different seasons through the spectral mean method; the statistical features include the standard deviation B1 of the B1-B7 bands in the original time series spectrum std 、B2 std 、B3 std 、B4 std 、B5 std 、B6 std 、B7 std , Standard deviation of Normalized Difference Vegetation Index NDVI std , Standard deviation of modified normalized water index MNDWI std , Standard deviation of normalized building index NDBI std ; The auxiliary feature is the annual mean value of night lights NTL Annual , the annual average value of night lights NTL Annual It is the annual night light data generated by the monthly night light data through the annual average synthesis algorithm.

[0026] The index features are generated by multi-source remote sensing data through the time series feature fusion method, including the normalized difference spectral index NDSI and the modified impervious surface index MISI. The calculation formula is as follows:

[0027] NDSI=ρm -ρ n

[0028] ρ m +ρ n

[0029] MISI = (1 - NDVI) × Log(NTL Annual + 1)

[0030] Among them, ρ m and ρ n are the spectral values of the m - band and n - band of the MODIS data respectively, where m, n = 1,..., 7, m ≠ n, and NTL Annual is the annual average value of the night - time light NTL Annual .

[0031] Preferably, in step S4, an impervious surface coverage estimation model is constructed. The impervious surface coverage estimation model is trained using the time - series impervious surface coverage samples and the corresponding impervious surface coverage feature sets. The impervious surface coverage feature sets of each year in the monitoring area are input into the impervious surface coverage estimation model, and the impervious surface coverage of the study area for each year is output.

[0032] Preferably, in step S5, the Mann - Kendall (MK) trend test is performed pixel - by - pixel on the impervious surface coverage results of each pixel during the monitoring time. For each time - series pixel, if P < 0.05 and Z ≥ 2.58, the change trend of the impervious surface coverage of this time - series pixel is increasing; if P < 0.05 and Z ≤ - 2.58, the change trend of the impervious surface coverage of this time - series pixel is decreasing; if P < 0.05 and - 2.58 < Z < 2.58, the change trend of the impervious surface coverage of this time - series pixel is unchanged; if P ≥ 0.05, the change trend of the impervious surface coverage of the time - series pixel is other; time - series correction is carried out based on the change situation and the impervious surface coverage results of each year.

[0033] The time - series correction uses a preset - length time window to judge the pixel trend. For the time - series pixels that conform to the change trend, the results of each year are retained. For the time - series pixels that do not conform to the change trend, time - series optimization is carried out.

[0034] The specific method of time - series optimization is as follows:

[0035] For the time - series pixels with an increasing trend that do not conform to the change trend, that is, the impervious surface coverage ISP t in the current year is greater than the impervious surface coverage ISP t+1 in the next year, calculate the average impervious surface coverage ISP mean within the preset - length time window. If ISP mean < ISP t+1, then ISP mean Assign value to ISP t If ISP mean >ISP t+1 , then ISP t+1 Assign value to ISP t ;

[0036] For the time series pixels that do not conform to the change trend in the reduction trend, that is, the impervious surface coverage ISP of the current year t Less than the impervious surface coverage ISP in the next year t+1 , calculate the mean impervious surface coverage ISP within a preset length time window mean If ISP mean >ISP t+1 , then ISP mean Assign value to ISP t If ISP mean <ISP t+1 , then ISP t+1 Assign value to ISP t ;

[0037] For time series pixels with unchanged trend test results, the average of all years’ impervious surface coverage is taken and assigned to the impervious surface coverage of each year;

[0038] For time series pixels with trend test results other than ISP, optimize year by year from later years to earlier years and calculate the mean impervious surface coverage ISP within the specified preset length time window. mean , the ISP mean Assign a value to the current year's impervious surface coverage ISP t .

[0039] The advantages of the present invention compared with the prior art are:

[0040] The present invention provides a large-scale dynamic impervious surface coverage estimation method based on time series optimization, which can be used to extract time-series impervious surface information in a large coverage area. The method obtains multi-temporal remote sensing data for monitoring impervious surface changes through the GEE platform, making full use of the advantages of MODIS data such as large coverage, long monitoring period and free acquisition, and giving full play to the ability of the GEE platform to store and process massive data, thereby improving the feasibility of large-scale impervious surface monitoring.

[0041] The present invention provides a large-scale dynamic impervious surface coverage estimation method based on time series optimization, which can be used to extract time series impervious surface information in a large coverage area. The method produces multiple types of time series impervious surface coverage samples based on multi-source remote sensing data, and constructs an impervious surface coverage remote sensing feature set. The features constructed by multi-source feature fusion and time series feature fusion technology include seasonal features, statistical features, auxiliary features, and index features. The features constructed by the present invention significantly reduce the confusion between impervious surfaces and other surface coverage, and the fused impervious surface coverage feature set has a positive effect on extracting impervious surface information in large-scale complex areas.

[0042] The present invention provides a large-scale dynamic impervious surface coverage estimation method based on time series optimization, which can be used to extract time series impervious surface information in a large coverage area. The method provides a time series impervious surface correction strategy to cope with complex changing trends, and corrects the impervious surface information of each year according to the increase, decrease, unchanged, and other trends of the time series impervious surface. This correction strategy improves the rationality of the extraction of time series impervious surface information. BRIEF DESCRIPTION OF THE DRAWINGS

[0043] The above characteristics, features and advantages of the present invention and their implementation methods and methods will become more clearly understood in conjunction with the following description of the embodiments, which are described in detail in conjunction with the accompanying drawings. Herein, a schematic diagram is shown:

[0044] Figure 1 The present invention provides an embodiment of a large-scale dynamic impervious surface coverage estimation method flow chart based on time series optimization;

[0045] Figure 2 Provide a flow chart of stable impervious surface coverage based on product derived timing in an embodiment of the present invention;

[0046] Figure 3 Provide a flowchart of timing optimization and correction in an embodiment of the present invention;

[0047] Figure 4 A schematic diagram of four types of time-series impervious surface change trends in the monitoring area in an embodiment of the present invention is provided;

[0048] Figure 5 This is a comparison diagram before and after pixel scale optimization in an embodiment of the present invention;

[0049] Figure 6 The present invention provides a schematic diagram of the distribution of impervious surface coverage in the monitoring area from 2000 to 2023 in the embodiment;

[0050] Figure 7 A schematic diagram of the distribution of impervious surface coverage in the local monitoring area from 2000 to 2023 is provided for the present invention in the embodiment;

[0051] Figure 8 A graph of the accuracy changes in the monitoring area from 2000 to 2023 in an embodiment of the present invention is provided. DETAILED DESCRIPTION

[0052] The present invention will be further described in detail below in conjunction with the accompanying drawings and embodiments, but the present invention is not limited in any way. Any changes or improvements made based on the present invention fall within the protection scope of the present invention.

[0053] This application proposes a large-scale dynamic impervious surface coverage estimation method based on time series optimization. Figure 1 , including the following steps:

[0054] Specifically, in this embodiment, taking the Yangtze River Delta and the surrounding areas of the North China Plain in China as examples, long-term temporal dynamic impervious surface coverage information from 2000 to 2023 is extracted.

[0055] S1: Acquisition and preprocessing of multi-source remote sensing data in the monitoring area: Through the GEE platform, long-term Landsat and MODIS data in the monitoring area are obtained, and outliers are removed using quality assessment bands to obtain time-series night light NTL data and multi-source impervious surface product data.

[0056] Remote sensing images of the monitoring area in each year are collected through the GEE platform, including Landsat surface reflectance (Landsat5 / 7 / 8) and MODIS surface reflectance (MOD09A1), forming a long-term remote sensing image sequence. Outliers are detected pixel by pixel based on the quality assessment bands of Landsat and MODIS, and abnormal pixels are masked.

[0057] The GEE platform is used to obtain Landsat5 / 7 / 8 time-series satellite remote sensing images including TM, ETM+, and OLI sensors with a spatial resolution of 30m. The Landsat5 / 7 / 8 time-series remote sensing dataset from 2000 to 2023 is selected within the monitoring area.

[0058] The GEE platform is used to obtain the MOD09A1 time series satellite remote sensing images of surface reflectance, with a spatial resolution of 500m and 46 periods per year. The MOD09A1 time series remote sensing dataset from 2000 to 2023 is selected within the monitoring area. MOD09A1 is an 8-day synthetic surface reflectance product of MODIS. Pixels are selected based on the high observation range, low viewing angle, no cloud or cloud shadow, and aerosol within 8 days. Each surface reflectance pixel contains the best Level-2G observation value. This product covers MODIS bands 1-7 as well as MODIS data quality layers and observation bands. Therefore, MOD09A1 has a large coverage area and can monitor a lot of time, which can ensure the detection of large-scale time series dynamic impervious surface coverage.

[0059] The multi-source remote sensing data used must be preprocessed before use to ensure that the coordinate system and spatial resolution of the remote sensing data source remain consistent.

[0060] According to the quality information recorded pixel by pixel in the StateQA file, including cloud layer, cloud shadow, rain and snow, the pixels that need to be removed and retained are marked, and the affected pixels are removed through the mask method to obtain high-quality time series observation data.

[0061] Multi-source impervious surface products contain products with different spatial and temporal resolutions, which are converted into impervious surface coverage with consistent temporal span and spatial resolution.

[0062] Obtain impervious surface data through the GEE platform and Zenodo platform, including global artificial impervious surface annual dynamic data products (GAIA, Global Artificial Impervious Areas), global impervious surface products (GISA, global impervious surface area), global impervious surface dynamic dataset (GISD, global impervious-surface dynamic dataset), global annual urban dynamics data (GAUD, global annualurban dynamics), global human settlement layer products (GHSL, Global Human Settlement Layer).

[0063] In terms of time series, multi-source remote sensing impervious surface products are annual and five-year interval data. The land cover types of adjacent years are similar. For the five-year interval data, the adjacent year filling method is used to supplement the blank time data until all existing data from 2000 to 2023 are covered.

[0064] A local grid with a spatial resolution of 500 m is set, and each grid represents the coverage within the range. For different products in each year, the area of ​​impervious surface within the 500 m grid is calculated to generate the impervious surface coverage with a spatial resolution of 500 m.

[0065] The multi-source impervious surface products are converted into impervious surface coverage that is consistent with the time span and spatial resolution of the remote sensing image. The calculation formula for the impervious surface coverage C is as follows:

[0066]

[0067] Among them, sum(AS) is the area of ​​the impervious surface in the grid, and length is the length of the grid, that is, the length of the pixel side.

[0068] S2: Production of time series impervious surface coverage samples: Obtain time series impervious surface coverage samples based on two methods, including the export of impervious surface products and time series sample supplementation, including the following steps:

[0069] S21: Obtain temporally stable impervious surface coverage samples based on impervious surface products. Obtain temporally stable impervious surface coverage samples through screening methods based on multi-source impervious surface products. Calculate the temporal fluctuation range of the same product at the pixel scale and the difference in impervious surface coverage between products. And identify stable pixels as temporally stable impervious surface coverage samples through the threshold method, such as Figure 2 shown.

[0070] Determine the sample screening area. Use the 500m impervious surface coverage converted from the five products of the same year to form the impervious surface coverage dataset for each year. Obtain the maximum coverage on the same pixel through the maximum value synthesis method based on the pixel scale. If the maximum coverage is greater than 0, there is an impervious surface; if the maximum coverage is equal to 0, there is no impervious surface, and these areas are eliminated using masks.

[0071] Sample screening. The existing impervious surface data has high reliability, but there are errors in the region. Therefore, obtaining samples with high similarity based on multi-source products is an effective way to improve sample accuracy. For each year, there are 5 products of impervious surface coverage on a pixel. The difference between the 5 products is calculated, and the product with a value less than the threshold is defined as a reliable product. The threshold is 0.05.

[0072] D ij =|AS i -AS j |

[0073] Among them, D ij It is the difference in impervious surface coverage between impervious surface product i and impervious surface product j.

[0074] At the pixel scale, the average value of impervious surface coverage of the screened multi-source products is calculated as a regional stable sample.

[0075] Extraction of time-series stable samples. Time-series stable samples refer to areas where the surface does not change over a long time series. Further screening of regional stable samples is performed. For regionally stable pixels, the coverage on the time series is obtained. The maximum and minimum values ​​of the impervious surface coverage from 2000 to 2023 are calculated on the pixel scale, and pixels with a range less than the threshold are selected, and the time series average value of the pixel is calculated as the time series stable sample. The threshold is set to 0.05.

[0076] S22: Supplement the temporal impervious surface coverage samples through manual correction, including temporal impervious surface coverage samples and non-impervious surface coverage samples in the same area.

[0077] The diversity and balance of input samples affect the prediction effect of the model, and the uniformity of spatial distribution and the diversity of sample types increase the prediction accuracy of the model. The above steps obtain samples with high reliability and similarity based on multi-source remote sensing data, but there are still deficiencies in the remaining types of samples.

[0078] Manual drawing supplements are made by adding time series samples based on medium and high-resolution remote sensing images and historical images provided by Google and Landsat.

[0079] S3: Construction of impervious surface coverage feature set: Based on multi-source remote sensing data, the impervious surface coverage feature set is constructed through multi-source feature fusion and time series feature fusion.

[0080] Based on the time series multi-source remote sensing data MOD09A1 and NTL, a feature set for estimating impervious surface coverage is constructed, including seasonal features, statistical features, auxiliary features, and index features. The feature types are shown in Table 1.

[0081] Table 1 Feature types

[0082]

[0083] Seasonal features are synthesized by time series MOD09A1 through time series fusion method. B1, B2, B3, B4, B5, B6, and B7 are synthesized by averaging algorithm from spectral bands of 46 MOD09A1 every year, including red band B1, near infrared band B2, blue band B3, green band B4, near infrared band B5, near infrared band B6, and short-wave infrared band B7. MOD09A1 images in summer and winter are selected respectively to generate spectral time series synthesis features of different seasons through spectral averaging method.

[0084] The estimation of impervious surface coverage usually selects summer as the information extraction season, because the vegetation in this season is significant and can better distinguish and identify impervious surfaces and vegetation. In large-scale information extraction, the seasonal changes of ground objects are diverse and complex. Therefore, obtaining summer and winter features separately can improve the richness of the feature set.

[0085] Statistical features include the standard deviation of the B1-B7 band in the time series original spectrum std 、B2 std 、B3 std 、B4 std 、B5 std 、B6 std 、B7 std , Standard deviation of Normalized Difference Vegetation Index (NDVI) std , Standard deviation of modified normalized water index MNDWI std , Standard deviation of normalized building index NDBI std . The Normalized Difference Vegetation Index is an index used to identify surface vegetation, which can increase the distinction between impervious surfaces and vegetation. The Modified Normalized Water Index is a commonly used index for identifying surface water bodies, which can increase the distinction between impervious surfaces and water bodies. The Normalized Building Index is used to identify buildings and built-up areas in remote sensing images, and has a positive effect on the extraction of impervious surface information. The more obvious the Normalized Building Index is, the greater the coverage of impervious surfaces.

[0086]

[0087] Among them, ρ red , green , swir , nir These are the red, blue, shortwave infrared, and near infrared bands in the MOD09A1 image.

[0088] The auxiliary feature is the annual mean value of night lights NTL Annual , the annual average value of night lights NTL Annual It is the annual night light data generated by the monthly night light data through the annual average synthesis algorithm. Night light data is closely related to human production activities and has potential value in large-scale impervious surface applications. Therefore, night light features make up for the limitations of spectral features in expressing impervious surface information.

[0089] The index features include the Normalized Difference Spectral Index (NDSI) and the Modified Impervious Surface Index (MISI). The Normalized Difference Spectral Index is synthesized by normalizing seven spectral features in pairs, which can fully express optical differences. Spectral information reflects the basic characteristics of surface changes in remote sensing images, and limited spectral information cannot meet the complex and changeable changes in surface information. However, pairwise normalization of spectral features can obtain more spectral differences, and the potential value of features can be mined through machine learning or deep learning methods. The Modified Impervious Surface Index is a remote sensing index synthesized by the Normalized Vegetation Index and night light features for extracting large-scale impervious surfaces. It comprehensively considers the spatial distribution of vegetation and night lights and is used to estimate the coverage of impervious surfaces. NDVI is negatively correlated with impervious surfaces, and NTL is positively correlated with impervious surfaces.

[0090] The normalized difference spectral index NDSI and the modified impervious surface index MISI are calculated by the following formulas:

[0091]

[0092] MISI=(1-NDVI)×Log(NTL Annual +1)

[0093] Among them, ρ m and ρ n are the spectral values ​​of the m-band and n-band of MODIS data, m,n=1,...,7,m≠n, NTL Annual is the annual average value of night lights NTL Annual .

[0094] S4: Construction of impervious surface coverage estimation model and extraction of spatiotemporal information: The temporal impervious surface coverage samples and the corresponding impervious surface coverage feature set are input into the impervious surface coverage estimation model for training, and the temporal impervious surface coverage results are estimated.

[0095] In this step, the impervious surface coverage estimation model is constructed using random forest regression, support vector machine regression, and neural network regression methods respectively. The time series impervious surface coverage samples and the impervious surface coverage feature set corresponding to the samples are used to train the impervious surface coverage estimation model. The impervious surface coverage feature set of the monitoring area for each year is input into the impervious surface coverage estimation model, and the impervious surface coverage of the study area for each year is output, that is, the large-scale spatiotemporal distribution results of impervious surface coverage are obtained for each year, and the value of each pixel represents the proportion of impervious surface in each area.

[0096] The annual impervious surface coverage is extracted using a regression algorithm, and the samples and features of each year are input into the prediction model. There are differences between different models, and the input samples and features are different. The feature scheme or model with the highest accuracy is selected as the annual impervious surface coverage result.

[0097] S5: Temporal logic optimization of impervious surface coverage: Perform temporal change correction based on the temporal impervious surface coverage under four impervious surface coverage change trends. The four impervious surface coverage change trends include: increase, decrease, unchanged and others.

[0098] The temporal logic optimization of impervious surface coverage is to correct the temporal impervious surface coverage according to different situations. The temporal impervious surface coverage changes under different situations refer to the long-term unchanged, increased, decreased and other changes of the impervious surface at the pixel scale. The temporal correction is carried out based on the changes and the annual impervious surface coverage results.

[0099] The main manifestation of the urbanization process is the process of permeable surface transformation into impermeable surface. At the same time, it is also the main form of human transformation of the earth's surface. Therefore, the time series prediction results of impermeable surface coverage can be optimized based on this long-term trend change. The impermeable surface coverage results are generated by independent models year by year, and there will be certain errors in the prediction results of each pixel year by year. Therefore, it is necessary to carry out time series optimization of the year-by-year results. The previous method formulated optimization rules based on the assumption that the impermeable surface is increasing or unchanged. However, the impermeable surface may increase or decrease in local areas, which is inevitable in time changes.

[0100] The MK trend test is performed pixel by pixel on the impervious surface coverage results of each pixel during the monitoring period to determine the overall change trend of the time series impervious surface coverage. There are four types of impervious surface coverage change trends: increase, decrease, unchanged, and other. The increase and decrease types represent that the impervious surface of the pixel increases or decreases as a whole during the monitoring period. The unchanged type represents that the impervious surface does not change during the monitoring period. The other types represent that the impervious surface changes unstably and fluctuates during the monitoring period.

[0101] In the MK test, the P value and Z value are indicators used to determine whether the data has a significant trend. The P value is the probability value obtained in the test, which represents the probability of the data showing the current trend when there is no trend. The Z value is the standardized value of the test statistic, which represents the strength of the data trend, as shown in Table 2.

[0102] Table 2 Temporal trends of impervious surface coverage

[0103]

[0104] According to the trend type, the time series results of impervious surface coverage are optimized year by year using a time window of length 3. For time series pixels that meet the trend of change, the results of each year are retained; for time series pixels that do not meet the trend of change, time series optimization is carried out. The order of time series optimization starts from the last year and gradually optimizes to the previous years. The results of the later years are more reliable.

[0105] The results of the independent prediction model are optimized in time series. The optimization steps are as follows: Figure 3 The specific steps are as follows:

[0106] (1) Increasing trend:

[0107] The results of the current year are in line with the trend of change and do not need to be revised, that is, the current year's impervious surface coverage ISP t Less than the impervious surface coverage ISP in the next year t+1 .

[0108] The results of the current year do not conform to the trend of change and need to be revised, that is, the current year's impervious surface coverage ISP t Greater than the impervious surface coverage ISP in the next year t+1 Calculate the mean impervious surface coverage (ISP) based on the time window mean If the ISP mean Smaller than ISP t+1 , ISP mean Alternative ISPs t ; If ISP mean Bigger than ISP t+1 , ISP t+1 Alternative ISPs t .

[0109]

[0110] Where N is the time window length, which is equal to 3, ISP mean is the average impervious surface cover calculated over the time window.

[0111] (2) Decreasing trend:

[0112] The results of the current year are in line with the trend of change and do not need to be revised, that is, the current year's impervious surface coverage ISP t Greater than the impervious surface coverage ISP in the next year t+1 .

[0113] The results of the current year do not conform to the trend of change and need to be revised, that is, the current year's impervious surface coverage ISP t Less than the impervious surface coverage ISP in the next year t+1Calculate the mean impervious surface coverage ISP based on the time window mean If the ISP mean Bigger than ISP t+1 , ISP mean Alternative ISPs t ; If ISP mean Smaller than ISP t+1 , ISP t+1 Alternative ISPs t .

[0114] (3) Unchanged trend:

[0115] The trend test result is unchanged, and the average impervious surface coverage ISP is synthesized using the results of all years. mean , replacing each year's results with the composite results.

[0116] (4) Other trends:

[0117] The mean impervious surface coverage ISP is synthesized using a time window of length 3 mean , use ISP mean Replace the result ISP of the current year t , optimizing the impervious surface coverage year by year from later years to earlier years.

[0118] The change trend of impervious surface coverage in the monitoring area is as follows: Figure 4 Most urban areas show an overall increasing trend, and the overall impervious surface coverage develops from low coverage to high coverage.

[0119] The optimization results of impervious surface coverage at the pixel scale in the monitoring area are as follows: Figure 5 As shown in the figure. The impervious surface in this area has reversed, the buildings have become vegetation, and the impervious surface has become permeable. According to the trend test results, the pixel trend change type is judged to be a reduction type, so the reduction type optimization method is used to carry out correction. After correction, the results of the next year are always equal to or less than the results of the previous year.

[0120] The time series of impervious surface coverage in the monitoring area is as follows: Figure 6 As shown in the figure, the impervious surface coverage of the monitoring area showed an overall expansion trend from 2000 to 2023. Figure 7 As shown in Figure 2, the inner part of the city tends to be stable, while the expansion of impervious surfaces occurs in the outer and scattered areas of the city.

[0121] S6: Result verification: Verify the accuracy of the extraction results based on the accuracy assessment method.

[0122] The accuracy of the obtained impervious surface coverage correction results was verified, and the goodness of fit R was selected. 2The accuracy of the extracted results is evaluated by the root mean square error (RMSE). Figure 8 shown.

[0123] The accuracy of the estimation of impervious surface coverage from 2000 to 2023 shows an upward trend. The recent data is of high quality and the impervious surface is widely distributed, so the accuracy is higher than in previous years. Among the three models, the random forest model has the best effect, the highest overall goodness of fit, and the smallest overall root mean square error. The support vector machine regression model has the worst prediction results.

[0124] In summary, large-scale impervious surface changes are a long-term process, and it is necessary to obtain long-term impervious surface coverage for detailed research on environmental monitoring and urban development. The present invention provides a large-scale dynamic impervious surface coverage estimation method based on time series optimization, which better extracts the evolution process of impervious surfaces under long-term changes in the past. The present invention constructs an impervious surface coverage feature set based on multi-source remote sensing data, inputs time series impervious surface samples into a regression prediction model, and optimizes the impervious surface coverage results of each year according to the overall change trend of the impervious surface. This method can improve the extraction accuracy and time series rationality of large-scale time series impervious surface coverage.

[0125] The present invention is not limited to the precise structure described above and shown in the drawings, and various modifications and changes can be made without departing from the scope thereof. Any modification, equivalent substitution, improvement, etc. made within the spirit and principle of the present invention shall be included in the protection scope of the present invention.

Claims

1. A large-scale dynamic impervious surface coverage estimation method based on time series optimization, characterized in that: The following steps are involved: S1: Acquisition and preprocessing of multi-source remote sensing data in the monitoring area: obtain long-term Landsat data and MODIS data in the monitoring area and remove outliers, obtain time-series night light NTL data and multi-source impervious surface product data; S2: Production of time series impervious surface coverage samples: Obtain time series impervious surface coverage samples through two methods: exporting impervious surface products and supplementing time series samples; S3: Construction of impervious surface coverage feature set: Based on the time series MODIS data and time series nighttime light NTL data, the impervious surface coverage estimation feature set is constructed, including seasonal features, statistical features, auxiliary features and index features; S4: Construction of impervious surface coverage estimation model and extraction of spatiotemporal information: Input the temporal impervious surface coverage samples and the corresponding impervious surface coverage feature set into the impervious surface coverage estimation model for training, and estimate the temporal impervious surface coverage results; S5: Temporal logic optimization of impervious surface coverage: carry out temporal change correction according to the temporal impervious surface coverage under four impervious surface coverage change trends, the four impervious surface coverage change trends include: increase, decrease, unchanged and others; S6: Result verification: Verify the accuracy of the extraction results according to the accuracy assessment method; In step S3, the seasonal features include summer features and winter features, which are synthesized by the time series MOD09A1 through the time series fusion method, and specifically include the red band B1, the near infrared band B2, the blue band B3, the green band B4, the near infrared band B5, the near infrared band B6, and the short-wave infrared band B7. The MOD09A1 images of summer and winter are respectively selected to generate the spectral time series synthesis features of different seasons through the spectral mean method; Statistical features include the standard deviation of the B1-B7 band in the time series original spectrum std 、B2 std 、B3 std 、B4 std 、B5 std 、B6 std 、B7 std , Standard deviation of Normalized Difference Vegetation Index NDVI std , Standard deviation of modified normalized water index MNDWI std , Standard deviation of normalized building index NDBI std ; The auxiliary feature is the annual mean value of night lights NTL Annual , the annual average value of night lights NTL Annual It is the annual night light data generated by the monthly night light data through the annual mean synthesis algorithm; The index features are generated by multi-source remote sensing data through a time series feature fusion method, including the normalized difference spectral index NDSI and the modified impervious surface index MISI.

2. According to the large-scale dynamic impervious surface coverage estimation method based on time series optimization according to claim 1, it is characterized in that: In step S1, remote sensing images of the monitoring area in each year are collected through the GEE platform, including Landsat surface reflectance data Landsat5 / 7 / 8 and MODIS surface reflectance data MOD09A1, to form a long-term remote sensing image sequence, and outliers are detected pixel by pixel according to the quality assessment bands of Landsat and MODIS, and abnormal pixels are masked; The multi-source impervious surface products are obtained through the GEE platform and the Zenodo platform, including: global artificial impervious surface annual dynamic data products, global impervious surface products, global impervious surface dynamic data sets, global annual urban dynamic data, and global human settlement layer products. The multi-source impervious surface products are converted into impervious surface coverage consistent with the time span and spatial resolution of remote sensing images. The calculation formula for impervious surface coverage C is as follows: Among them, sum(AS) is the area of ​​the impervious surface in the grid, and length is the length of the grid, that is, the length of the pixel side.

3. According to the large-scale dynamic impervious surface coverage estimation method based on time series optimization according to claim 1, it is characterized in that: In step S2, the preparation of the time series impervious surface coverage sample includes the following steps: S21: Extract time-stable impervious surface coverage samples from the impervious surface product; S22: Supplement the time-series impervious surface coverage samples through manual correction, including the impervious surface coverage samples with time-series changes in the same area and the non-impervious surface coverage samples.

4. The large-scale dynamic impervious surface coverage estimation method based on time series optimization according to claim 3 is characterized in that: In step S21, time-series stable impervious surface coverage samples are obtained from multi-source impervious surface products through the threshold method. If the time-series fluctuation range of the impervious surface coverage of the same product and the difference in impervious surface coverage between different products are within the preset threshold range at the pixel scale, they are identified as stable pixels and used as time-series stable impervious surface coverage samples.

5. The large-scale dynamic impervious surface coverage estimation method based on time series optimization according to claim 1 is characterized in that: The calculation formulas of the Normalized Difference Spectral Index (NDSI) and the Modified Impervious Surface Index (MISI) are as follows: DAY=(1-NDVI)×Log(NTL Annual +1) Among them, ρ m and ρ n are the spectral values ​​of the m-band and n-band of MODIS data, m,n=1,...,7,m≠n, NTL Annual is the annual average value of night lights NTL Annual .

6. The large-scale dynamic impervious surface coverage estimation method based on time series optimization according to claim 1 is characterized in that: In step S4, an impervious surface coverage estimation model is constructed. The time-series impervious surface coverage samples and the corresponding impervious surface coverage feature sets of the samples are used to train the impervious surface coverage estimation model. The impervious surface coverage feature sets of each year in the monitoring area are input into the impervious surface coverage estimation model, and the impervious surface coverage of each year in the study area is output.

7. The large-scale dynamic impervious surface coverage estimation method based on time series optimization according to claim 1 is characterized in that: In step S5, the Mann-Kendall (MK) trend test is performed pixel by pixel on the impervious surface coverage results of each pixel during the monitoring time. For each time-series pixel, if P < 0.05 and Z ≥ 2.58, the change trend of the impervious surface coverage of this time-series pixel is increasing; if P < 0.05 and Z ≤ -2.58, the change trend of the impervious surface coverage of this time-series pixel is decreasing; if P < 0.05 and -2.58 < Z < 2.58, the change trend of the impervious surface coverage of this time-series pixel remains unchanged; if P ≥ 0.05, the change trend of the impervious surface coverage of the time-series pixel is other; time-series correction is carried out based on the change situation and the impervious surface coverage results of each year.

8. The large-scale dynamic impervious surface coverage estimation method based on time series optimization according to claim 7 is characterized in that: Time-series correction uses a preset-length time window to judge the pixel trend. For time-series pixels that conform to the change trend, the results of each year are retained. For time-series pixels that do not conform to the change trend, time-series optimization is carried out.

9. The large-scale dynamic impervious surface coverage estimation method based on time series optimization according to claim 8 is characterized in that: The specific method of time-series optimization is as follows: For the time series pixels that do not conform to the changing trend in the increasing trend, that is, the impervious surface coverage ISP of the current year t Greater than the impervious surface coverage ISP in the next year t+1 , calculate the mean impervious surface coverage ISP within a preset length time window mean If ISP mean <ISP t+1 , then ISP mean Assign value to ISP t If ISP mean >ISP t+1 , then ISP t+1 Assign value to ISP t ; For the time series pixels that do not conform to the change trend in the reduction trend, that is, the impervious surface coverage ISP of the current year t Less than the impervious surface coverage ISP in the next year t+1 , calculate the mean impervious surface coverage ISP within a preset length time window mean If ISP mean >ISP t+1 , then ISP mean Assign value to ISP t If ISP mean <ISP t+1 , then ISP t+1 Assign value to ISP t ; For time-series pixels with an unchanged trend test result, the average value of the impervious surface coverage of all years is taken and assigned to the impervious surface coverage of each year; For time series pixels with trend test results other than ISP, optimize year by year from later years to earlier years and calculate the mean impervious surface coverage ISP within the specified preset length time window. mean , the ISP mean Assign a value to the current year's impervious surface coverage ISP t .

Images