CN110929739B - Automatic impervious surface range remote sensing iterative extraction method - Google Patents

Abstract

The invention discloses an automatic remote sensing iterative extraction method of impermeable surface range, which realizes the automatic extraction of impermeable surface coverage information in a continuous large area by integrating nighttime light images and Landsat TM images. The present invention firstly extracts urban areas, surrounding areas and non-urban areas with different aggregation densities of impermeable surfaces from nighttime light images, and then performs spectral index statistics and threshold segmentation corresponding to different areas on the Landsat TM image to obtain impervious The initial classification samples of surfaces/non-impervious surfaces are used to train the classifier. Finally, in the classification process, iteratively integrates spectral information and spatial information to automatically select new classification samples until the classification results are stable, and the impervious surface range information extraction process is completed. The invention can realize automatic extraction of impermeable surface range information with stable precision and high efficiency.

Description

An Automatic Remote Sensing Iterative Extraction Method of Impervious Surface Range

technical field

The invention belongs to the technical field of remote sensing information extraction, and in particular relates to the design of an automatic remote sensing iterative extraction method for an impermeable surface range.

Background technique

Impervious surfaces are special types of land cover that prevent water and air from penetrating into the ground. They are usually made of man-made materials, including roads, building roofs, and parking lots. The expansion of impervious surfaces can be regarded as an intuitive indicator of the degree of urbanization development, and it also has a significant impact on the urban regional environment. For example, the impact on watershed non-point source pollution; the impact on urban heat island effect, etc. Therefore, accurate and reliable impervious distribution information and its time-series change information play an important and fundamental role in urban research and management. The extraction and monitoring of impervious surfaces is an important and active research direction in the field of remote sensing.

At present, a variety of methods based on remote sensing have been widely used in the extraction of impervious surface information. The main research methods are: (1) Mixed pixel decomposition, according to the "V-I-S" model, the heterogeneous urban land cover is simplified to vegetation, The impervious surface and the soil are composed of three types of surface features, and the impervious surface can be defined as an independent end member or as the sum of high reflectivity end members and low reflectivity end members, so that according to the least square method and other mixed (2) Statistical regression model, by classifying the high-resolution remote sensing images to obtain the real impervious surface reference data on the ground, and then establishing a statistical regression model between it and the characteristic bands of the multispectral remote sensing images Solve the proportion of impervious surface in the pixel; (3) Spectral index, through the combination of remote sensing image bands, construct the spectral index that can strengthen the information of impervious surface and weaken other information, and then obtain the impervious surface area through threshold segmentation; (4 ) Image classification model, using conventional classifiers such as classification decision trees, artificial neural networks, support vector machines, or deep learning models such as convolutional neural networks, to classify remote sensing images and extract impervious surface information.

In addition, the analysis method using multi-source remote sensing data and multi-type feature integration has also been applied to the research of impervious surface extraction. For example, the National Land Cover Database (NLCD) of the United States uses thematic data to assist in the extraction of impervious surface products.

To sum up, most of the existing research needs to manually obtain a certain number of impervious surface samples to train or adjust the parameters of the extraction model. Due to the variety of impervious surfaces and their spectral complexity, and other ground objects (such as dry soil, exposed rocks, etc.) have serious spectral confusion, it is difficult to automate the extraction of large-scale research through effective spectral indices The extraction process, while the manual sampling method is inefficient and cannot guarantee the representativeness of the sample.

Contents of the invention

The purpose of the present invention is to realize spatially continuous extraction of impermeable surface information in a larger area, and to solve the deficiency of relying on artificial samples in existing methods, and propose an automatic remote sensing iterative extraction method of impermeable surface range.

The technical solution of the present invention is: an automatic remote sensing iterative extraction method for impermeable surface range, comprising the following steps:

S1. Obtain DSMP/OLS nighttime light remote sensing data and Landsat TM5 multispectral remote sensing data of the area to be extracted, and obtain DMSP/OLS images and TM images of the area to be extracted.

S2. According to the DMSP/OLS image of the region to be extracted, the city region is extracted to obtain the city region, the surrounding region of the city and the non-urban region.

S3. According to the TM image of the area to be extracted, samples are automatically selected for the urban area and the non-urban area respectively, and the ISA sample set of the urban area and the non-ISA sample set of the non-urban area are obtained to form an initial classification sample set.

S4. Use a classifier to iteratively classify the initial classification sample set, and automatically collect new classification samples by integrating spectral information and spatial information during the iterative process, and extract ISAs of urban areas, urban peripheral areas, and non-urban areas.

S5. Integrate the ISAs of urban areas, urban peripheral areas and non-urban areas, and complete the ISA extraction of the area to be extracted.

Further, step S2 includes the following sub-steps:

S21. According to the DMSP/OLS image of the area to be extracted, obtain the light area whose night light brightness value is equal to 0 in the area to be extracted, and use it as a non-urban area.

S22. According to the DMSP/OLS image of the area to be extracted, obtain the light area in the area to be extracted whose nighttime light brightness value is greater than 0, and convert gas into light isosurface objects of different brightness according to the light brightness value.

S23. Arrange all light isosurface objects in descending order according to their area size.

S24. Starting from the light isosurface object with the largest area, sequentially detect whether all light isosurface objects contain other light isosurface objects, and if so, proceed to step S25; otherwise, the light area corresponding to the light isosurface object Set to city area.

S25. Set the light area corresponding to the light iso-surface object that only has a containment relationship among the included light iso-surface objects as an urban area.

S26. Remove the urban area from the light area to obtain the surrounding area of the city.

Further, step S3 includes the following sub-steps:

S31. Calculate the NDVI index, NDWI index and NDBI index of each pixel on the TM image of the region to be extracted.

S32. Combine the NDVI index, NDWI index and NDBI index of each pixel to obtain a composite index of each pixel, and obtain a composite index image according to the composite index of each pixel.

S33. According to the eigenvalue and NDVI index of the pixel points in the non-urban area on the synthetic index image, the non-ISA samples of the non-urban area are obtained by threshold segmentation, and a preset number of non-ISA samples are selected to obtain an initial non-ISA sample set. .

S34. According to the feature value and NDVI index of the pixel points in the urban area on the synthetic index image, obtain the ISA samples of the urban area through threshold segmentation, and select a preset number of ISA samples to obtain an initial ISA sample set.

Further, the calculation formula of the composite index of each pixel in step S32 is:

Among them, Idx _i represents the composite index of the i-th pixel, and NDBI _i , NDVI _i , and NDWI _i represent the NDBI index, NDVI index and NDWI index of the i-th pixel, respectively.

Further, the segmentation formulas of non-ISA samples in non-urban areas in step S33 and the segmentation formulas of ISA samples in urban areas in step S34 are:

Among them, Pixel _i represents the i-th pixel, ISA represents an ISA sample, soil represents a soil sample, vegetation represents a vegetation sample, and the soil sample and the vegetation sample together constitute a non-ISA sample, N/A represents an undetermined sample, and a _Idx represents a synthetic ISA segmentation threshold of the index, b _Idx represents the soil segmentation threshold of the synthetic index, c _NDVI represents the vegetation segmentation threshold of the NDVI index, U is the set of urban area pixels, and R is the set of non-urban area pixels.

Further, step S4 includes the following sub-steps:

S41. Input the initial classification sample set into the C5.0 decision tree classifier to classify the Landsat TM images in urban areas, and obtain the initial classification results and corresponding classification probabilities.

S42. On the initial classification result image, according to the classification probability thresholds of ISA and non-ISA types, classify the pixels smaller than the threshold as "unclassified" type, and retain the classification type corresponding to the pixels larger than the threshold, and obtain partial classification results image.

S43. On some of the classification result images, perform statistics on the "unclassified" pixels, and calculate their corresponding local spatial feature values.

S44. Add the local spatial feature value and classification probability value corresponding to the "unclassified" pixel by weight, sort the "unclassified" pixel from large to small according to the weighted addition result, and select the local area in each category The pixel with the highest spatial feature value is added to the classification sample set as a new sample.

S45. Determine whether the area difference between the ISA of the urban area in the current classification result and the ISA of the urban area in the previous iterative classification result is less than the preset error threshold, if so, obtain the ISA of the urban area, and enter step S46, otherwise return to step S45. S41 adopts a new sample set to iteratively classify.

S46. Taking the Landsat TM image of the surrounding area of the city as the input of the C5.0 decision tree classifier, performing iterative classification according to the same method as steps S41-S45, to obtain the ISA of the surrounding area of the city.

S47. Taking the Landsat TM image of the non-urban area as the input of the C5.0 decision tree classifier, performing iterative classification according to the same method as steps S41-S45, to obtain the ISA of the non-urban area.

Further, the calculation formula of the local spatial feature value in step S43 is:

Among them, Spa(w _k ) _i means that the i-th pixel corresponds to the local spatial feature value of the w _k-th class, and D _i,j means the Euclidean distance between the i-th pixel and the j-th class w _k classed pixel The distance, b, c are predefined distance contribution adjustment parameters.

Further, the calculation formula for the weighted addition of the local spatial feature value and the classification probability value in step S44 is:

Score(w _k ) _i ＝a ₁ ×Prob(w _k ) _i +a ₂ ×Spa(w _k ) _i

Among them, Score(w _k ) _i represents the weighted addition result of the local spatial feature value corresponding to the i-th pixel and the classification probability value, and Prob(w _k ) _i represents the classification probability of the i-th pixel corresponding to the w _k -th class value, Spa(w _k ) _i means that the i-th pixel corresponds to the local spatial feature value of the w _k -th class, and a ₁ and a ₂ are predefined weight factors.

The beneficial effects of the present invention are: the remote sensing iterative extraction method for impermeable surface range proposed by the present invention integrates the spatial and spectral information in night light images and Landsat TM images, and realizes the automatic extraction of ISA within the entire image range. It is an important improvement compared with the existing manual selection of ISA samples, and effectively improves the extraction efficiency of ISA. Experimental results prove that the present invention better solves the problem of spectral confusion between ISA and soil, and the extraction accuracy is high and stable; by comparing with NLCD test data, the average overall accuracy and kappa coefficient in urban areas are 88.23% and 88.23% respectively. 0.63; 78.6% and 0.54 outside the urban area, both of which are better than the extraction accuracy of the manual extraction method.

Description of drawings

FIG. 1 is a flow chart of an automatic remote sensing iterative extraction method for an impermeable surface range provided by an embodiment of the present invention.

Fig. 2 is a schematic diagram showing the location of the experimental area and the corresponding DSMP/OLS, Landsat TM image and verification area data provided by the embodiment of the present invention.

FIG. 3 is a schematic diagram of the results of extracting urban areas, surrounding areas and non-urban areas from DMSP/OLS images provided by an embodiment of the present invention.

Fig. 4 is a Landsat TM image map of a typical urban area provided by an embodiment of the present invention.

FIG. 5 is a schematic diagram of the extraction results of impermeable surface ranges in typical urban areas provided by the embodiment of the present invention.

Detailed ways

Exemplary embodiments of the present invention will now be described in detail with reference to the accompanying drawings. It should be understood that the implementations shown and described in the drawings are only exemplary, intended to explain the principle and spirit of the present invention, rather than limit the scope of the present invention.

An embodiment of the present invention provides a method for automatically extracting a large-area impermeable surface area (Impervious Surface Area, ISA), as shown in FIG. 1 , including the following steps S1-S5:

S1. Obtain DSMP/OLS (Defense Meteorological Satellite Program's Operational Line Scan System) night light remote sensing data and Landsat TM5 multispectral remote sensing data of the area to be extracted, and obtain the DMSP/OLS image and TM of the area to be extracted. image.

In the embodiment of the present invention, the row number of the TM image is Path15/Row30, and the image is composed of other 6 bands (visible light band, near-infrared band and two mid-infrared bands) except the thermal infrared band, because the atmospheric effect is very important for the Landsat image. Operations such as classification have little effect, so only the DN value of the image is corrected as the apparent reflectance of the atmosphere. For the night light data, the stable light data provided by DMSP/OLS is selected. This product removes the interference of flashes and clouds, and is often used to estimate the indicator factors of human activity areas in previous studies. Because the spatial resolution of the DMSP/OLS image is 30 arc seconds (about 1 km), in order to match the resolution of the TM image, it was reprojected to the UTM coordinate system and resampled to 30m and cropped to the same size as the TM image.

S2. According to the DMSP/OLS image of the region to be extracted, the city region is extracted to obtain the city region, the surrounding region of the city and the non-urban region.

On the DMSP/OLS image, the intensity of the light directly indicates the position information of the ISA accumulation area. However, urban areas of different scales have different intensity and range on night light images. In addition to the divergence of the light itself, the brightness of lights in different urban areas will also interfere with each other, so it is difficult to extract through a fixed threshold. In the embodiment of the present invention, urban areas are extracted by analyzing the spatial relationship between light areas with different intensities. First, in the light areas with luminance values greater than 0 in night light images, they are converted into isosurface objects with different luminance values according to the light luminance values. On this basis, urban areas are defined as isosurface objects that are spatially independent of each other, that is, they are required to meet the following conditions:

(1) The isosurface object does not contain other isosurface objects;

(2) If the isosurface object contains other isosurface objects, there is only a contained or contained spatial relationship between the contained objects, and there is no intersecting spatial relationship.

Step S2 includes the following sub-steps S21-S26:

S21. According to the DMSP/OLS image of the area to be extracted, obtain the light area whose night light brightness value is equal to 0 in the area to be extracted, and use it as a non-urban area;

S22. According to the DMSP/OLS image of the area to be extracted, obtain the light area in the area to be extracted whose night light brightness value is greater than 0, and convert gas into light isosurface objects of different brightness according to the light brightness value;

S23. Arranging all light isosurface objects in descending order according to their area size;

S24. Starting from the light isosurface object with the largest area, sequentially detect whether all light isosurface objects contain other light isosurface objects, and if so, proceed to step S25; otherwise, the light area corresponding to the light isosurface object set to city area;

S25. Set the light area corresponding to the light iso-surface object that only has a containment relationship among the included light iso-surface objects as an urban area;

S26. Remove the urban area from the light area to obtain the surrounding area of the city.

In the embodiment of the present invention, after processing the DMSP/OLS image, the entire image is divided into three parts: urban area, urban peripheral area and non-urban area.

S3. According to the TM image of the area to be extracted, samples are automatically selected for the urban area and the non-urban area respectively, and the ISA sample set of the urban area and the non-ISA sample set of the non-urban area are obtained to form an initial classification sample set.

Step S3 includes the following sub-steps S31-S34:

S31. Calculate the NDVI index (Normalized Difference Vegetation Index, normalized vegetation index), NDWI index (Normalized Difference Water Index, normalized water index) and NDBI index (Normalized Difference Building) of each pixel on the TM image of the area to be extracted Index, normalized building index).

S32. Combine the NDVI index, NDWI index and NDBI index of each pixel to obtain a composite index of each pixel, and obtain a composite index image according to the composite index of each pixel.

The formula for calculating the synthetic index of each pixel is:

Among them, Idx _i represents the composite index of the i-th pixel, and NDBI _i , NDVI _i , and NDWI _i represent the NDBI index, NDVI index and NDWI index of the i-th pixel, respectively.

S33. According to the eigenvalue and NDVI index of the pixel points in the non-urban area on the synthetic index image, the non-ISA samples of the non-urban area are obtained by threshold segmentation, and a preset number of non-ISA samples are selected to obtain an initial non-ISA sample set. .

S34. According to the feature value and NDVI index of the pixel points in the urban area on the synthetic index image, obtain the ISA samples of the urban area through threshold segmentation, and select a preset number of ISA samples to obtain an initial ISA sample set.

The segmentation formula of the non-ISA sample of the non-urban area in step S33 and the segmentation formula of the ISA sample of the urban area in the step S34 are:

Among them, Pixel _i represents the i-th pixel, ISA represents an ISA sample, soil represents a soil sample, vegetation represents a vegetation sample, and the soil sample and the vegetation sample together constitute a non-ISA sample, N/A represents an undetermined sample, and a _Idx represents a synthetic ISA segmentation threshold of the index, b _Idx represents the soil segmentation threshold of the synthetic index, c _NDVI represents the vegetation segmentation threshold of the NDVI index, U is the set of urban area pixels, and R is the set of non-urban area pixels.

S4. Use a classifier to iteratively classify the initial classification sample set, and automatically collect new classification samples by integrating spectral information and spatial information during the iterative process, and extract ISAs of urban areas, urban peripheral areas, and non-urban areas.

Step S4 includes the following sub-steps S41-S47:

S41. Input the initial classification sample set into the C5.0 decision tree classifier to classify the Landsat TM images in urban areas, and obtain the initial classification results and corresponding classification probabilities.

S42. On the initial classification result image, according to the classification probability thresholds of ISA and non-ISA types, classify the pixels smaller than the threshold as "unclassified" type, and retain the classification type corresponding to the pixels larger than the threshold, and obtain partial classification results image.

S43. On some classification result images, perform statistics on the "unclassified" pixels, and calculate their corresponding local spatial feature values, the calculation formula is:

Among them, Spa(w _k ) _i means that the i-th pixel corresponds to the local spatial feature value of the w _k-th class, and D _i,j means the Euclidean distance between the i-th pixel and the j-th class w _k classed pixel The distance, b, c are predefined distance contribution adjustment parameters.

S44. Add the local spatial feature value and classification probability value corresponding to the "unclassified" pixel by weight, sort the "unclassified" pixel from large to small according to the weighted addition result, and select the local area in each category The pixel with the highest spatial feature value is added to the classification sample set as a new sample.

The calculation formula for the weighted addition of local spatial feature values and classification probability values is:

Score(w _k ) _i ＝a ₁ ×Prob(w _k ) _i +a ₂ ×Spa(w _k ) _i

Among them, Score(w _k ) _i represents the weighted addition result of the local spatial feature value corresponding to the i-th pixel and the classification probability value, and Prob(w _k ) _i represents the classification probability of the i-th pixel corresponding to the w _k -th class value, Spa(w _k ) _i means that the i-th pixel corresponds to the local spatial feature value of the w _k -th class, and a ₁ and a ₂ are predefined weight factors.

S45. Determine whether the area difference between the ISA of the urban area in the current classification result and the ISA of the urban area in the previous iterative classification result is less than the preset error threshold, if so, obtain the ISA of the urban area, and enter step S46, otherwise return to step S45. S41 adopts a new sample set to iteratively classify.

S46. Taking the Landsat TM image of the surrounding area of the city as the input of the C5.0 decision tree classifier, performing iterative classification according to the same method as steps S41-S45, to obtain the ISA of the surrounding area of the city.

S47. Taking the Landsat TM image of the non-urban area as the input of the C5.0 decision tree classifier, performing iterative classification according to the same method as steps S41-S45, to obtain the ISA of the non-urban area.

S5. Integrate the ISAs of urban areas, urban peripheral areas, and non-urban areas to obtain the impermeable surface range extraction results of the entire area to be extracted, and complete the ISA extraction of the area to be extracted.

The effect of the remote sensing iterative extraction method for impermeable surface range provided by the present invention will be further described below with a set of specific experimental examples.

Experimental example 1:

Taking Chengdu, Sichuan Province, China as the area to be extracted, this area contains large-scale urban areas, covering the ISA surface from low density to high density, and is a representative area for impervious surface extraction.

The step S2 of the present invention is used to extract the urban area, and the urban area extraction mainly depends on the isosurface information of the light intensity in the DMSP/OLS image. In the experiment, the brightness interval is set to 10 to extract isovalue lines, and merged with the boundary line of the light area (threshold is 1), and these isovalue lines are converted into isosurfaces. After setting the minimum detection area to 9km ² (corresponding to the area of 30 TM image pixels), the results of urban area extraction in these isosurfaces are shown in Figure 3. A total of 41 urban areas were detected in the study area. Due to the low spatial resolution of the DMSP/OLS image itself and the divergent characteristics of the light, the boundaries of the detection results cannot completely match the boundaries of the cities on the TM image, which is more obvious than that on the TM image. The boundary of is larger, so it is difficult to evaluate the detection accuracy by quantitative comparison on the boundary morphology. Therefore, the detection result is compared with the corresponding online map and TM image by visual inspection. If the detection area contains city points on the map, the detection result is correct. Experimental results show that the detection accuracy of the urban area extraction algorithm provided by the present invention is greater than 95%, and can detect relatively small urban areas. After detection, the entire image is divided into three areas, namely: urban area, urban peripheral area and non-urban area.

Experimental example 2:

The step S3 of the present invention is used to make statistics on the spectral index images of each area and perform threshold segmentation, firstly divide the non-ISA type sample area in the non-urban area, and then divide the ISA type sample area in the urban area. In order to ensure high accuracy in sample selection, the corresponding band segmentation threshold is set as the sum of the average value and standard deviation of the index statistics in the area, and the thresholds such as a _Idx , b _Idx and c _NDVI are obtained for segmentation. In the outer urban areas, 500 samples of soil types and 500 samples of vegetation types were selected by the random sampling method in the corresponding segmentation results; Sampling method for sample selection.

Experimental example 3:

The samples obtained in the pre-order process were used as the initial training set and input into an iterative classification framework for training, in which the C5.0 decision tree classifier was used as the core classifier to classify the reflectance values of the 6 bands of the TM image into feature bands. In each iteration process, the classifier outputs the pixels whose classification probability is greater than 0.9 as the classified pixels of the corresponding category, and outputs other pixels as unclassified pixels; the local spatial feature calculation is performed on the unclassified pixels, and the The calculated value is normalized to [0, 1]; then the spatial feature value and classification probability are summed according to the weight of 0.5, and the 500 pixels with the largest value in each category are selected as new classification samples and added to the training set Carry out the next training until the number of iterations is greater than 15 or the difference between the area of the two ISA extraction results before and after is less than 10%, and then stop the iteration to obtain the ISA extraction result of this area. As a comparison, at the same time, manually select samples in this area, and then extract the ISA results through the decision tree classifier.

In the ISA extraction results of the test area, 15,000 test points are randomly selected and the corresponding test set is tested point-to-point. The extraction accuracy is shown in Table 1.

Table 1

It can be seen from Table 1 that the average overall accuracy and Kappa coefficient of the automatic remote sensing iterative extraction method for impermeable surface range provided by the present invention are 89.69% and 0.908, respectively, which are better than the accuracy of the manual sampling method. This is because most of the samples obtained by the manual sampling method are samples with clear spectral characteristics, and there is a lack of collection of samples of mixed pixel types, so it is easy to generate errors in the omission of ISA. The invention continuously selects new samples to add through an iterative process, so that the representativeness of the samples in terms of spatial distribution and spectral characteristics is better. Moreover, after the calculated value of the local spatial feature provided by the present invention is integrated into the new sample selection, the overall classification accuracy is improved, and the error standard deviation of the multiple experimental results is reduced. This shows that the present invention can provide statistically sound extraction and classification results, especially for the ISA producer accuracy and non-ISA user accuracy are significantly improved.

On this basis, three representative regions were selected to demonstrate the ISA extraction effect of the present invention. The TM images corresponding to these three areas are shown in Figure 4; the manual sampling classification results corresponding to each area are shown in Figure 5(a), the extraction results of 10 iterations are shown in Figure 5(b), and the extraction results of 15 iterations are shown in Figure 5( c). Compared with the TM image, it can be seen that the extraction result provided by the present invention shows more complete ISA distribution information and fewer ISA classification errors. This is because new samples generated during iterations increase the representativeness of the training set. In particular, comparing the two extraction results provided by the present invention, it can be seen that the result with more iterations has less noise information, especially the confusion between the better resolved ISA and the bare soil category, which is common in ISA extraction methods The problem. Therefore, the method provided by the present invention can realize an automatic ISA extraction process without manual intervention and reference data, and can provide higher accuracy than general manual sampling methods.

Those skilled in the art will appreciate that the embodiments described here are to help readers understand the principles of the present invention, and it should be understood that the protection scope of the present invention is not limited to such specific statements and embodiments. Those skilled in the art can make various other specific modifications and combinations based on the technical revelations disclosed in the present invention without departing from the essence of the present invention, and these modifications and combinations are still within the protection scope of the present invention.

Claims (7)

1. an automatic impermeable surface range remote sensing iterative extraction method, is characterized in that, comprises the following steps: S1. Obtain DSMP/OLS nighttime light remote sensing data and Landsat TM5 multispectral remote sensing data of the area to be extracted, and obtain DMSP/OLS images and TM images of the area to be extracted; S2. According to the DMSP/OLS image of the area to be extracted, the urban area is extracted to obtain the urban area, the surrounding area of the city and the non-urban area; S3. According to the TM image of the area to be extracted, samples are automatically selected for the urban area and the non-urban area respectively, and the ISA sample set of the urban area and the non-ISA sample set of the non-urban area are obtained to form an initial classification sample set; S4. Use a classifier to iteratively classify the initial classification sample set, and automatically collect new classification samples by integrating spectral information and spatial information during the iterative process, and extract the ISAs of urban areas, urban peripheral areas and non-urban areas; S5. Integrating the ISAs of the urban area, the surrounding area of the city and the non-urban area, and completing the extraction of the ISA of the area to be extracted; Described step S4 comprises following sub-steps: S41. Input the initial classification sample set into the C5.0 decision tree classifier to classify the Landsat TM image in the urban area, and obtain the initial classification result and the corresponding classification probability; S42. On the initial classification result image, according to the classification probability thresholds of ISA and non-ISA types, classify the pixels smaller than the threshold as "unclassified" type, and retain the classification type corresponding to the pixels larger than the threshold, and obtain partial classification results image; S43. On some of the classification result images, perform statistics on the "unclassified" pixels, and calculate their corresponding local spatial feature values; S44. Add the local spatial feature value and classification probability value corresponding to the "unclassified" pixel by weight, sort the "unclassified" pixel from large to small according to the weighted addition result, and select the local area in each category The pixel with the highest spatial feature value is added to the classification sample set as a new sample; S45. Determine whether the area difference between the ISA of the urban area in the current classification result and the ISA of the urban area in the previous iterative classification result is less than the preset error threshold, if so, obtain the ISA of the urban area, and enter step S46, otherwise return to step S45. S41 adopting a new sample set for iterative classification; S46. Using the Landsat TM image of the surrounding area of the city as the input of the C5.0 decision tree classifier, perform iterative classification according to the same method as steps S41 to S45, to obtain the ISA of the surrounding area of the city; S47. Taking the Landsat TM image of the non-urban area as the input of the C5.0 decision tree classifier, performing iterative classification according to the same method as steps S41-S45, to obtain the ISA of the non-urban area. 2. The impermeable surface range remote sensing iterative extraction method according to claim 1, wherein said step S2 comprises the following sub-steps: S21. According to the DMSP/OLS image of the area to be extracted, obtain the light area whose night light brightness value is equal to 0 in the area to be extracted, and use it as a non-urban area; S22. According to the DMSP/OLS image of the area to be extracted, obtain the light area in the area to be extracted whose night light brightness value is greater than 0, and convert gas into light isosurface objects of different brightness according to the light brightness value; S23. Arranging all light isosurface objects in descending order according to their area size; S24. Starting from the light isosurface object with the largest area, sequentially detect whether all light isosurface objects contain other light isosurface objects, and if so, proceed to step S25; otherwise, the light area corresponding to the light isosurface object set to city area; S25. Set the light area corresponding to the light iso-surface object that only has a containment relationship among the included light iso-surface objects as an urban area; S26. Remove the urban area from the light area to obtain the surrounding area of the city. 3. The impermeable surface range remote sensing iterative extraction method according to claim 1, wherein said step S3 comprises the following sub-steps: S31. Calculate the NDVI index, NDWI index and NDBI index of each pixel on the TM image of the area to be extracted; S32. Combining the NDVI index, NDWI index and NDBI index of each pixel to obtain a composite index of each pixel, and obtain a composite index image according to the composite index of each pixel; S33. According to the eigenvalue and NDVI index of the pixel points in the non-urban area on the synthetic index image, the non-ISA samples of the non-urban area are obtained by threshold segmentation, and a preset number of non-ISA samples are selected to obtain an initial non-ISA sample set. ; S34. According to the feature value and NDVI index of the pixel points in the urban area on the synthetic index image, obtain the ISA samples of the urban area through threshold segmentation, and select a preset number of ISA samples to obtain an initial ISA sample set. 4. the impermeable surface range remote sensing iterative extraction method according to claim 3, is characterized in that, the computing formula of the synthetic index of each pixel in the described step S32 is:

Among them, Idx _i represents the composite index of the i-th pixel, and NDBI _i , NDVI _i , and NDWI _i represent the NDBI index, NDVI index and NDWI index of the i-th pixel, respectively. 5. the impermeable surface range remote sensing iterative extraction method according to claim 4, is characterized in that, the segmentation formula of the non-ISA sample of the non-urban area in the described step S33 and the segmentation formula of the ISA sample of the urban area in the step S34 are both for:

Among them, Pixel _i represents the i-th pixel, ISA represents an ISA sample, soil represents a soil sample, vegetation represents a vegetation sample, and the soil sample and the vegetation sample together constitute a non-ISA sample, N/A represents an undetermined sample, and a _Idx represents a synthetic ISA segmentation threshold of the index, b _Idx represents the soil segmentation threshold of the synthetic index, c _NDVI represents the vegetation segmentation threshold of the NDVI index, U is the set of urban area pixels, and R is the set of non-urban area pixels. 6. the impermeable surface range remote sensing iterative extraction method according to claim 1, is characterized in that, the computing formula of local space characteristic value in the described step S43 is:

Among them, Spa(w _k ) _i means that the i-th pixel corresponds to the local spatial feature value of the w _k-th class, and D _i,j means the Euclidean distance between the i-th pixel and the j-th class w _k classed pixel The distance, b, c are predefined distance contribution adjustment parameters. 7. The impermeable surface range remote sensing iterative extraction method according to claim 1, characterized in that, the calculation formula of the weighted addition of the local spatial feature value and the classification probability value in the step S44 is: Score(w _k ) _i ＝a ₁ ×Prob(w _k ) _i +a ₂ ×Spa(w _k ) _i Among them, Score(w _k ) _i represents the weighted addition result of the local spatial feature value corresponding to the i-th pixel and the classification probability value, and Prob(w _k ) _i represents the classification probability of the i-th pixel corresponding to the w _k -th class value, Spa(w _k ) _i means that the i-th pixel corresponds to the local spatial feature value of the w _k -th class, and a ₁ and a ₂ are predefined weight factors.

Images