A remote sensing detection method for urban construction land change

Abstract

The invention discloses a method for remote sensing detection of urban construction land change, comprising: inputting two Landsat5-TM remote sensing images of different phases that have been geometrically corrected and radiometrically corrected; performing geometric precise registration on the two remote sensing images; performing principal component analysis dimensionality reduction processing on the seven bands of the two remote sensing images, and obtaining the corresponding first components respectively; calculating the texture features of the TM7 band of the two remote sensing images using a gray level co-occurrence matrix; calculating the normalized vegetation index NDVI of the two remote sensing images; calculating the normalized building index NDBI of the two remote sensing images; performing image difference calculation on the obtained four pairs of change detection factors; performing principal component analysis on the obtained four difference result gray images, and taking the first group thereof as the urban construction land change detection result; using a threshold to segment the urban construction land change detection result gray image, and extracting a binary image for change detection. The invention has the technical advantages of simple and easy algorithm, stable and reliable, efficient calculation, and accurate calculation.

Description

Technical Field

[0001] The present invention relates to the technical field of remote sensing satellite image processing, and in particular to a land satellite-oriented urban construction land change remote sensing detection calculation method, device and electronic equipment. Background Art

[0002] Land cover change is the result of the combined effects of human activities and natural processes on the earth's land surface, and is an important part of global change research. Cities are areas with the most active human activities and the highest degree of land use. The transformation of land use types from natural types such as cultivated land, grassland, and forest land to artificial construction land is the most important form of transformation of urban land use during urbanization. Urban expansion is the spread of urban areas in space. It is characterized by the continuous expansion of construction land. It occurs widely and rapidly around the world, and has a profound impact on resources, environment, geographical processes, biodiversity, and other aspects. Remote sensing detection of changes in urban construction land can provide a scientific basis for urban management, urban planning, and ecological and environmental research, provide data support for the scientific development of cities, and has important practical significance for the construction of ecological cities.

[0003] The Landsat satellite series is jointly managed by NASA and the United States Geological Survey (USGS). Since 1972, the LANDSAT satellite series has been launched one after another. It is a series of earth observation satellite systems used by the United States to detect earth resources and environment. The LANDSAT satellite series is equipped with the Thematic Mapper (TM) or Enhanced Thematic Mapper (ETM+), which passively senses the solar radiation reflected by the surface and the thermal radiation emitted. However, there is no method for calculating urban construction land based on Landsat satellite remote sensing that is stable, reliable, efficient and accurate through simple and easy algorithms. Summary of the invention

[0004] The technical problem to be solved by the present invention is to solve the deficiencies of the above-mentioned prior art and to provide a method and an implementation device for remote sensing detection of urban construction land changes based on Landsat satellites, which is accurate in calculation, simple and easy to operate, stable and reliable, and efficient in operation.

[0005] The technical solution adopted by the present invention to solve its technical problem is:

[0006] A calculation method for remote sensing detection of urban construction land change includes the following steps:

[0007] S1: Input two Landsat5-TM remote sensing images of different phases, which have been geometrically corrected and radiometrically corrected, and obtain seven bands from TM1 to TM7;

[0008] S2: Use remote sensing image A as a reference to perform geometric precise registration of remote sensing image B;

[0009] S3: Perform principal component analysis and dimensionality reduction processing on seven bands from TM1 to TM7 of remote sensing images A and B respectively, and obtain seven components from PC1 to PC7. The first component (PC1) corresponding to remote sensing images A and B is taken as the change detection factors F1A and F1B respectively;

[0010] S4: Perform texture analysis on TM7 of remote sensing images A and B respectively to obtain the corresponding texture feature mean statistic grayscale image, ie, change detection factors F2A and F2B;

[0011] S5: Calculate the normalized vegetation index of remote sensing images A and B respectively to obtain their corresponding normalized vegetation index, i.e., change detection factors F3A and F3B;

[0012] S6: Calculate the normalized building indexes of the remote sensing images A and B respectively, and obtain their corresponding normalized building indexes, i.e., change detection factors F4A and F4B;

[0013] S7: Based on the obtained four pairs of change detection factors F1A and F1B, F2A and F2B, F3A and F3B, and F4A and F4B, perform image difference operation on each pair of factors to obtain corresponding difference operation results DF1, DF2, DF3, and DF4;

[0014] S8: The four difference result grayscale images of DF1, DF2, DF3 and DF4 are combined into a new four-band image C by overlapping the bands, and the four bands of image C are subjected to principal component analysis to obtain a total of four components from PC1 to PC4, and the first component (PC1) is taken as the grayscale image D of the urban construction land change detection result;

[0015] S9: For the grayscale image D of the urban construction land change detection result, the change detection threshold T is obtained using the image histogram, and the pixels with grayscale values ​​greater than the threshold T are set as changed pixels through band operation to obtain the urban construction land change binary image E.

[0016] As a preferred solution of the present invention, in S1: the remote sensing image selects a level 4 product, and all seven bands thereof are resampled to a spatial resolution of 30 meters.

[0017] As a preferred solution of the present invention, in S2, the root mean square error of the registration is required to be controlled within 0.5 pixels.

[0018] As a preferred solution of the present invention, in S4, the texture feature mean statistic grayscale image is obtained by calculating the grayscale co-occurrence matrix based on a 3×3 window size, wherein the grayscale co-occurrence matrix is ​​defined as: assuming that a certain area of ​​the image has N grayscale values, the grayscale co-occurrence matrix corresponding to the area is an N×N order probability matrix, and the element at position (i,j) (1, …, i, … N; 1, …, j, …, N) in the matrix is ​​the probability P(i,j) of the phenomenon that the pixel with grayscale i leaves a certain fixed position relationship δ=(DX, DY) and the pixel with grayscale j is j, δ is called the displacement, DX and DY are called the distance;

[0019] The feature quantities used for texture recognition are some statistics calculated from the gray-level co-occurrence matrix. Based on the gray-level co-occurrence matrix, the calculation formula for the mean statistics of texture features is:

[0020] ;

[0021] Among them, P(i,j) is the value of the gray-level co-occurrence matrix, and i and j represent the pixel coordinates of a point in the image space.

[0022] As a preferred embodiment of the present invention, the calculation formula of the normalized vegetation index is: NDVI = (TM4-TM3) / (TM4+TM3), wherein NDVI is the normalized vegetation index, TM3 is the grayscale value of the third band of the Landsat5-TM remote sensing image, and TM4 is the grayscale value of the fourth band of the Landsat5-TM remote sensing image.

[0023] As a preferred solution of the present invention, the calculation formula of the normalized building index is: NDBI= (TM5-TM4) / (TM5+ TM4), wherein, NDBI is the normalized building index, TM4 is the grayscale value of the 4th band of the Landsat5-TM remote sensing image, and TM5 is the grayscale value of the 5th band of the Landsat5-TM remote sensing image.

[0024] A device for implementing a calculation method for remote sensing detection of urban construction land changes includes the following modules:

[0025] (1) Image acquisition module, used to input at least two Landsat5-TM remote sensing images of different phases that have been geometrically corrected and radiometrically corrected;

[0026] (2) Image precision registration module, used to perform geometric precision registration of remote sensing image B using remote sensing image A as a reference;

[0027] (3) Image dimension reduction processing module, used to perform principal component analysis on the seven bands of remote sensing images A and B to obtain the corresponding first components respectively;

[0028] (4) Image texture analysis module, used to calculate the texture features of the TM7 band of remote sensing images A and B using the gray-level co-occurrence matrix

[0029] (5) Vegetation index calculation module, used to calculate the normalized difference vegetation index NDVI of remote sensing images A and B respectively;

[0030] (6) Building index calculation module, used to calculate the normalized building index NDBI of remote sensing images A and B respectively;

[0031] (7) An image difference operation module, used to perform image difference operation on four pairs of change detection factors;

[0032] (8) A principal component analysis module is used to combine the four difference result grayscale images obtained into a new four-band image by overlapping the bands, and perform principal component analysis on the image, taking the first group as the urban construction land change detection result;

[0033] (9) Threshold segmentation module, which is used to segment the grayscale image of the urban construction land change detection result using the threshold and extract the corresponding binary image.

[0034] The beneficial effects of the present invention are as follows: the urban construction land change detection calculation method and device provided by the present invention have the calculation advantages of simple and easy algorithm, stable and reliable, efficient operation and accurate calculation. BRIEF DESCRIPTION OF THE DRAWINGS

[0035] Figure 1 It is a flow chart of a method for remote sensing detection and calculation of urban construction land change in an embodiment of the present invention;

[0036] Figure 2 is a grayscale image of the change detection factor F1 in an embodiment of the present invention;

[0037] Figure 3 is a grayscale image of the change detection factor F2 in an embodiment of the present invention;

[0038] Figure 4 is a grayscale image of the change detection factor F3 in an embodiment of the present invention;

[0039] Figure 5 is a grayscale image of the change detection factor F4 in an embodiment of the present invention;

[0040] Figure 6 Grayscale images of the difference calculation results of the change detection factors DF1(a), DF2(b), DF3(c), and DF4(d) in the embodiment of the present invention;

[0041] Figure 7 The grayscale image D(a) of the urban construction land change detection result and the binary image E(b) after threshold separation in the embodiment of the present invention;

[0042] Figure 8 It is a schematic diagram of the structure of a device for implementing the remote sensing detection and calculation method of urban construction land change in an embodiment of the present invention. DETAILED DESCRIPTION

[0043] The present invention is further described below in conjunction with the accompanying drawings and examples, but these specific implementation schemes do not limit the protection scope of the present invention in any way.

[0044] Example 1

[0045] refer to Figure 1-Figure 7 :

[0046] S1, remote sensing image input

[0047] Two Landsat5-TM remote sensing images A and B at different time phases that have been geometrically corrected and radiometrically corrected are input. Furthermore, the two Landsat5-TM remote sensing images are selected as Level 4 products, and all seven bands thereof are resampled to a spatial resolution of 30 meters, and the corresponding first to seventh bands are TM1, TM2, TM3, TM4, TM5, TM6, and TM7, respectively. Specifically, the embodiment of the present invention selects Level 4 (Level 4 products) Landsat5-TM remote sensing images of Changsha area at two time phases, namely, October 4, 1996 (A) and September 6, 2009 (B), that have been geometrically corrected and radiometrically corrected, and cuts out a part of the area as the input data of the embodiment of the present invention;

[0048] S2, precise registration of remote sensing images

[0049] The remote sensing image A on October 4, 1996 is used as the reference to perform geometric precise registration of the remote sensing image B on September 6, 2009. The root mean square error of the registration is required to be controlled within 0.5 pixels.

[0050] S3, dimensionality reduction processing of remote sensing images

[0051] The principal component analysis was performed on the 7 bands TM1 to TM7 of remote sensing images A and B respectively, and 7 components PC1 to PC7 were obtained. The first component (PC1) corresponding to remote sensing images A and B was taken as the change detection factors F1A and F1B respectively. Specifically, Figure 2 a. Figure 2 b, respectively, are the change detection factors F1A and F1B calculated by the land satellite in the embodiment of the present invention;

[0052] S4, Remote Sensing Image Texture Analysis

[0053] The grayscale co-occurrence matrix is ​​calculated based on a 3×3 window size to obtain the grayscale image of the mean statistic of the texture feature. The grayscale co-occurrence matrix is ​​defined as: Assuming that there are N grayscale values ​​in a certain area of ​​the image, the grayscale co-occurrence matrix corresponding to the area is an N×N probability matrix. The element at position (i,j) (1, …, i, … N; 1, …, j, …, N) in the matrix is ​​the probability P(i,j) that the pixel with grayscale i leaves a fixed position relationship δ=(DX, DY) and the pixel with grayscale j is j. δ is called the displacement, DX and DY are called distances, and the feature quantities for texture recognition are some statistics calculated from the grayscale co-occurrence matrix. Based on the grayscale co-occurrence matrix, the calculation formula for the mean statistic of the texture feature is:

[0054]

[0055] Among them, P(i,j) is the value of the gray-level co-occurrence matrix, i and j represent the pixel coordinates of a point in the image space, and then the TM7 of the remote sensing images A and B are texture analyzed to obtain the corresponding texture feature mean statistic grayscale image, that is, the change detection factors F2A and F2B. Specifically, Figure 3 a. Figure 3 b, respectively, are the change detection factors F2A and F2B calculated by the land satellite in the embodiment of the present invention;

[0056] S5, Normalized Difference Vegetation Index calculation

[0057] The calculation formula of normalized difference vegetation index is:

[0058] NDVI= (TM4- TM3) / (TM4+ TM3),

[0059] Among them, NDVI is the Normalized Difference Vegetation Index,

[0060] TM3—the grayscale value of the third band of the Landsat5-TM remote sensing image.

[0061] TM4—Grayscale value of the 4th band of Landsat5-TM remote sensing image,

[0062] Calculate the normalized vegetation index of remote sensing images A and B respectively to obtain their corresponding normalized vegetation index, that is, change detection factors F3A and F3B. Specifically, Figure 4 a. Figure 4 b, respectively, are the change detection factors F3A and F3B calculated by the land satellite in the embodiment of the present invention;

[0063] S6, Normalized building calculation

[0064] The calculation formula of the normalized building index is:

[0065] NDBI= (TM5- TM4) / (TM5+ TM4),

[0066] Among them, NDBI—Normalized Difference Building Index,

[0067] TM4—Grayscale value of the 4th band of Landsat5-TM remote sensing image,

[0068] TM5—Grayscale value of the 5th band of Landsat5-TM remote sensing image,

[0069] Calculate the normalized building index of remote sensing images A and B respectively, and obtain their corresponding normalized building index, that is, change detection factors F4A and F4B. Specifically, Figure 5 a. Figure 5 As shown in b, they are respectively the change detection factors F4A and F4B calculated by the land satellite in the embodiment of the present invention.

[0070] S7, image difference calculation

[0071] Based on the obtained 4 pairs of change detection factors F1A and F1B, F2A and F2B, F3A and F3B, F4A and F4B, each pair of factors is subjected to image difference operation to obtain corresponding difference operation results DF1, DF2, DF3 and DF4. Specifically, Figure 6 a. Figure 6 b. Figure 6 c and Figure 6 d, which are grayscale images DF1, DF2, DF3 and DF4 of the difference calculation results of the change detection factors calculated by the land satellite in the embodiment of the present invention;

[0072] S8, principal component analysis

[0073] The four difference grayscale images of DF1, DF2, DF3 and DF4 are combined into a new four-band image C by overlapping the bands. The four bands of image C are subjected to principal component analysis to obtain four components from PC1 to PC4. The first component (PC1) is taken as the grayscale image D of the urban construction land change detection result. Specifically, Figure 7 As shown in a, it is a grayscale image D of the urban construction land change detection result calculated by the land satellite in an embodiment of the present invention;

[0074] S9, threshold segmentation

[0075] For the grayscale image D of the urban construction land change detection result, the threshold T of change detection is obtained using the image histogram. The pixels with grayscale values ​​greater than the threshold T are set as changed pixels through band operation to obtain the urban construction land change binary image E. Specifically, Figure 7 As shown in b, it is a binary image E after the urban construction land change threshold is separated by land satellite calculation in an embodiment of the present invention.

[0076] The detection result of the change of urban construction land in Changsha calculated by the embodiment of the invention is consistent with the actual situation, and the conclusion is:

[0077] Due to urban expansion in Changsha County, large tracts of agricultural land and grassland around the central urban area have been converted into urban construction land; large tracts of wetlands on the north side of Binhe South Road have been converted into residential land, with major residential areas including Vanke City, Zhujiang Licheng and Wanguo City; the city in the Hexi Lu Valley-Meixi Lake area has expanded rapidly westward; large tracts of urban construction land have been added around Xiangjiang North Road and Changsha North Station due to the construction of Jinxia Cross-border Industrial and Trade City.

[0078] More specifically, a method similar to classification accuracy evaluation is used to evaluate the accuracy of the results of urban construction land change detection. First, the change detection result is used as a classification image to include two categories: "change" and "no change"; then, by selecting test samples of the "change" category and the "no change" category on the original remote sensing image, the accuracy evaluation is achieved by establishing a confusion matrix. In the embodiment of the invention, the overall detection accuracy of the urban construction land change detection result is 91.28%, the Kappa coefficient is 0.6917, the false detection rate is 10.95%, and the missed detection rate is 12.87%. The urban construction land change has achieved a high detection accuracy, showing the superiority of the invention detection method. Therefore, the urban construction land change detection based on this method has excellent calculation accuracy and strong stability.

[0079] It can be seen from the above description that the urban construction land change remote sensing detection calculation method and electronic equipment for land satellites of the present invention have technical advantages such as simple and easy algorithm, stable and reliable, high operation efficiency and accurate calculation. From the perspective of market application, it is easy to form functional modules and embed them into various remote sensing image processing software, and has broad application prospects and market potential.

[0080] On the other hand, the embodiment of the present invention also provides a device for implementing the remote sensing detection method for urban construction land changes. Figure 8 , is a schematic diagram of the structure of a device for implementing the urban construction land change remote sensing detection method according to an embodiment of the present invention. It includes:

[0081] The image acquisition module 101 inputs two Landsat5-TM remote sensing images of different phases that have been geometrically corrected and radiometrically corrected;

[0082] An image precise registration module 102 performs geometric precise registration of remote sensing image B based on remote sensing image A;

[0083] The image dimension reduction processing module 103 performs principal component analysis on the seven bands of the remote sensing images A and B to obtain the corresponding first components respectively;

[0084] The image texture analysis module 104 calculates the texture features of the TM7 band of the remote sensing images A and B using the gray level co-occurrence matrix;

[0085] The vegetation index calculation module 105 calculates the normalized vegetation index NDVI of the remote sensing images A and B respectively;

[0086] A building index calculation module 106 calculates the normalized building index NDBI of the remote sensing images A and B respectively;

[0087] An image difference operation module 107 performs image difference operation on four pairs of change detection factors;

[0088] The principal component analysis module 108 combines the four obtained difference result grayscale images into a new four-band image by overlapping the bands, and performs principal component analysis, and takes the first group thereof as the urban construction land change detection result;

[0089] The threshold segmentation module 109 segments the grayscale image of the urban construction land change detection result using the threshold value to extract the corresponding binary image.

[0090] The above-mentioned embodiment is a remote sensing detection and calculation device for urban construction land changes, which is used to implement the corresponding remote sensing detection and calculation method for urban construction land changes in the above-mentioned embodiment, and has the beneficial effects of the corresponding method embodiment, which will not be repeated here.

[0091] It should be understood by those skilled in the art that the discussion of any of the above embodiments is merely illustrative and is not intended to imply that the scope of the present disclosure (including the claims) is limited to these examples. Within the scope of the present invention, the technical features in the above embodiments or different embodiments may also be combined, and there are many other variations of different aspects of the present invention as above, which are not provided in detail for the sake of simplicity.

[0092] The embodiments of the present invention are intended to cover all such substitutions, modifications and variations that fall within the broad scope of the appended claims. Therefore, any omissions, modifications, equivalent substitutions, improvements, etc. made within the spirit and principles of the present invention should be included in the protection scope of the present invention.

Claims

1. A calculation method for remote sensing detection of urban construction land changes, characterized in that: It includes the following steps: S1: Input two Landsat5-TM remote sensing images of different phases, which have been geometrically corrected and radiometrically corrected, and obtain seven bands from TM1 to TM7; S2: using one of the two remote sensing images as a reference to perform geometric precise registration on the other remote sensing image; S3: performing principal component analysis dimensionality reduction processing on the seven bands of the two remote sensing images respectively, outputting seven components, and taking the first group corresponding to the remote sensing image as the first pair of change detection factors; S4: performing texture analysis on TM7 of the two remote sensing images respectively, and obtaining the corresponding texture feature mean statistic grayscale images as the second pair of change detection factors; S5: respectively calculating the normalized vegetation indexes of the two remote sensing images, and obtaining the corresponding normalized vegetation indexes as the third pair of change detection factors; S6: respectively calculating the normalized building indexes of the two remote sensing images, and obtaining the corresponding normalized building indexes as a fourth pair of change detection factors; S7: Based on the four pairs of change detection factors obtained in S3-S6, perform image difference operation on each pair of factors to obtain corresponding four pairs of difference operation results; S8: The four grayscale images of the difference operation results obtained in S7 are combined into a new four-band image by overlapping the bands, and the four bands of the image are subjected to principal component analysis to output four components, and the first component is taken as the grayscale image of the urban construction land change detection result; S9: For the grayscale image of the urban construction land change detection result described in S8, the threshold of change detection is obtained using the image histogram, and the pixels with grayscale values ​​greater than the threshold are set as changed pixels through band operation to obtain a binary image of urban construction land change.

2. The calculation method for remote sensing detection of urban construction land changes according to claim 1 is characterized in that: In S1, the remote sensing image selects a level 4 product, and all seven bands thereof are resampled to a spatial resolution of 30 meters.

3. The calculation method for remote sensing detection of urban construction land changes according to claim 1 is characterized by: In S2, the root mean square error of the registration is required to be controlled within 0.5 pixels.

4. The calculation method for remote sensing detection of urban construction land changes according to claim 1 is characterized in that: In S4, the texture feature mean statistic grayscale image is obtained by calculating the grayscale co-occurrence matrix based on a 3×3 window size. The grayscale co-occurrence matrix is ​​defined as follows: assuming that a certain area of ​​the image has N grayscale values, the grayscale co-occurrence matrix corresponding to the area is an N×N probability matrix, and the element at position (i,j) (1, …, i, … N; 1, …, j, …, N) in the matrix is ​​the probability P(i,j) of the phenomenon that the pixel with grayscale i leaves a certain fixed position relationship δ=(DX, DY) and the pixel with grayscale j is located there, δ is called the displacement, and DX and DY are called the distance; The feature quantities used for texture recognition are some statistics calculated from the gray-level co-occurrence matrix. Based on the gray-level co-occurrence matrix, the calculation formula for the mean statistics of texture features is: ； Among them, P(i,j) is the value of the gray-level co-occurrence matrix, and i and j represent the pixel coordinates of a point in the image space.

5. The calculation method for remote sensing detection of urban construction land changes according to claim 1 is characterized in that: The calculation formula of the normalized vegetation index is: NDVI = (TM4-TM3) / (TM4+TM3), wherein NDVI is the normalized vegetation index, TM3 is the grayscale value of the third band of the Landsat5-TM remote sensing image, and TM4 is the grayscale value of the fourth band of the Landsat5-TM remote sensing image.

6. A calculation method for remote sensing detection of urban construction land changes according to any one of claims 1 to 5, characterized in that: The calculation formula of the normalized building index is: NDBI = (TM5-TM4) / (TM5+TM4), where NDBI is the normalized building index, TM4 is the grayscale value of the 4th band of the Landsat5-TM remote sensing image, and TM5 is the grayscale value of the 5th band of the Landsat5-TM remote sensing image.

7. A device for implementing the calculation method for remote sensing detection of urban construction land changes according to claim 1, characterized in that: It includes the following modules: (1) Image acquisition module, used to input at least two Landsat5-TM remote sensing images of different phases that have been geometrically corrected and radiometrically corrected; (2) Image precision registration module, used to perform geometric precision registration of remote sensing image B using remote sensing image A as a reference; (3) Image dimension reduction processing module, used to perform principal component analysis on the seven bands of remote sensing images A and B to obtain the corresponding first components respectively; (4) Image texture analysis module, used to calculate the texture features of TM7 band of remote sensing images A and B using gray-level co-occurrence matrix; (5) Vegetation index calculation module, used to calculate the normalized difference vegetation index NDVI of remote sensing images A and B respectively; (6) Building index calculation module, used to calculate the normalized building index NDBI of remote sensing images A and B respectively; (7) An image difference operation module, used to perform image difference operation on four pairs of change detection factors; (8) A principal component analysis module is used to combine the four difference result grayscale images obtained into a new four-band image by overlapping the bands, and perform principal component analysis on the image, taking the first group as the urban construction land change detection result; (9) Threshold segmentation module, which is used to segment the grayscale image of the urban construction land change detection result using the threshold and extract the corresponding binary image.

Images