Land utilization classification method based on space-time fusion technology

Abstract

The invention relates to the technical field of land classification, in particular to a land utilization classification method based on a space-time fusion technology, which comprises the following steps: firstly, acquiring original remote sensing images, classifying the original remote sensing images, and then respectively creating buffer areas for the divided original remote sensing images; the method comprises the following steps of: acquiring an original remote sensing image, dividing a buffer area into buffer area groups according to the category of the original remote sensing image, carrying out space-time fusion on different buffer area groups through a space-time fusion algorithm, splicing the fused image to generate a final image, and finally carrying out land utilization classification by adopting a supervised classification method according to the final image and a historical image. According to the method and the device, an effective solution is provided for improving the land utilization classification precision by integrating the time dynamic information of the low-resolution image and the space details of the high-resolution image based on the space-time fusion technology.

Description

Technical Field

[0001] The present invention relates to the technical field of land classification, and in particular to a land use classification method based on spatio-temporal fusion technology. Background Art

[0002] With the rapid development of remote sensing technology, land use classification has become a key technical means in the fields of resource management, environmental protection, agricultural monitoring, and urban planning. The spatial and temporal resolutions of remote sensing images have an important impact on the accuracy of land use classification. Low-resolution images have advantages in dynamic change monitoring due to their strong time series characteristics; high-resolution images are mainly characterized by clear spatial details and can better represent the boundaries of ground objects. However, relying solely on low-resolution or high-resolution images is difficult to balance temporal dynamics and spatial details, and it is difficult to meet the classification requirements of complex ground object areas (such as building boundaries, vegetation areas, etc.). Summary of the Invention

[0003] In view of this, the purpose of the present invention is to provide a land use classification method based on spatio-temporal fusion technology to overcome the problems existing in the current prior art.

[0004] To achieve the above purpose, the present invention adopts the following technical solutions:

[0005] The present application provides a land use classification method based on spatio-temporal fusion technology, including:

[0006] Obtain the original remote sensing image and classify the original remote sensing image;

[0007] Create buffer areas for the divided original remote sensing images respectively, and divide the buffer areas into buffer area groups according to the categories of the original remote sensing images;

[0008] Perform spatio-temporal fusion on different buffer area groups through a spatio-temporal fusion algorithm;

[0009] Stitch the fused images to generate the final image;

[0010] Perform land use classification using a supervised classification method according to the final image and historical images.

[0011] Further, in the above method, the obtaining the original remote sensing image and classifying the original remote sensing image includes:

[0012] Obtain the original remote sensing image;

[0013] Use the supervised classification method to compare the feature vectors of each pixel point P of the original remote sensing image with the training sample set to obtain the class label;

[0014] Classify the original remote sensing image according to the category label; wherein, the classification categories include: building area, vegetation area, and other areas;

[0015] Extract the building area and the vegetation area as vector files.

[0016] Further, for the method described above, create buffer areas for the divided original remote sensing image respectively, and divide the buffer areas into buffer area groups according to the categories of the original remote sensing image, including:

[0017] Dynamically adjust the buffer width of the building area according to the building density of the building area;

[0018] Dynamically adjust the buffer width of the vegetation area according to the vegetation coverage and time dynamics of the vegetation area;

[0019] According to the buffer width of the building area and the buffer width of the vegetation area, generate the buffer of the building area and the buffer of the vegetation area respectively through the buffer generation formula;

[0020] For the overlapping part of the building area buffer and the vegetation area buffer, retain the information of the high-priority area according to the preset priority, cut off the low-priority area, and introduce a fuzzy assignment mechanism at the junction of the building area buffer and the vegetation area buffer;

[0021] According to the cropped building area buffer and vegetation area buffer, use ArcGIS tools to create corresponding building buffer groups, corresponding vegetation buffer groups, and other area boundary groups, and convert them into SHP files.

[0022] Further, for the method described above, perform spatio-temporal fusion on different buffer area groups through a spatio-temporal fusion algorithm, including:

[0023] Use the STDF-A algorithm to perform spatio-temporal fusion on the building buffer group;

[0024] Use the FIT-FIC algorithm to perform spatio-temporal fusion on the vegetation buffer group;

[0025] Use the FSDAF algorithm to perform spatio-temporal fusion on the other area boundary group through a multi-scale feature extraction method.

[0026] Further, for the method described above, splice the fused images to generate the final image, including:

[0027] Determine the reference image and the adjustment image in the fused image;

[0028] Perform grayscale adjustment and stitching on the adjusted image according to the reference image;

[0029] Perform color normalization on the adjusted image after grayscale adjustment according to the reference image;

[0030] Perform boundary smoothing on the stitched fused image to generate the final image.

[0031] Further, for the method described above, the performing grayscale adjustment and stitching on the adjusted image according to the reference image includes:

[0032] Calculate the grayscale histograms of the reference image and the adjusted image respectively to determine the reference image grayscale histogram and the adjusted image grayscale histogram;

[0033] Calculate the cumulative distribution functions of the reference image grayscale histogram and the adjusted image grayscale histogram respectively;

[0034] Perform grayscale adjustment on the adjusted image according to the cumulative distribution function of the reference image grayscale histogram and the cumulative distribution function of the adjusted image grayscale histogram, and stitch the adjusted adjusted image and the reference image to complete the grayscale adjustment of the adjusted image.

[0035] Further, for the method described above, the performing color normalization on the adjusted image after grayscale adjustment according to the reference image includes:

[0036] Calculate the mean and standard deviation of the reference image and the adjusted image respectively;

[0037] Adjust the mean and standard deviation of the adjusted image according to the mean and standard deviation of the reference image through an adjustment formula to complete the color normalization of the adjusted image;

[0038] Adjust the brightness and contrast of the stitching area of the stitched image.

[0039] Further, for the method described above, the performing boundary smoothing on the stitched fused image to generate the final image includes:

[0040] Define the junction area for the stitched fused image;

[0041] Assign weights to the junction area;

[0042] Perform weighted summation of the pixel values of the junction area according to the assigned weights to complete the weighted average fusion of the fused image;

[0043] Construct a Gaussian pyramid for the image region after weighted average fusion, and decompose the image region into sub-images with multiple resolutions;

[0044] Generate a Laplacian pyramid in the Gaussian pyramid, and extract the detail information of sub-images with different resolutions;

[0045] Generate a weight pyramid corresponding to the Laplacian pyramid based on the weight calculation of the junction region;

[0046] Perform weighted fusion on the Gaussian pyramid, the Laplacian pyramid, and the weight pyramid at each resolution layer;

[0047] Reconstruct the fused pyramid from low resolution to high resolution in sequence to generate the final image.

[0048] Further, for the method described above, before performing land use classification using the supervised classification method based on the final image and historical images, it further includes:

[0049] Obtain a reference image;

[0050] Calculate the peak signal-to-noise ratio, structural similarity, and spectral consistency between the final image and the reference image respectively;

[0051] Evaluate the quality of the final image according to the peak signal-to-noise ratio, the structural similarity, and the spectral consistency.

[0052] Further, the method described above further includes:

[0053] Obtain the land use classification result and the ground truth data;

[0054] Calculate the overall classification accuracy of the land use classification result through the confusion matrix method;

[0055] Calculate the Kappa coefficient between the land use classification result and the ground truth data;

[0056] Calculate the user accuracy of the land use classification result and the producer accuracy of the ground truth data;

[0057] Evaluate the land use classification result according to the overall classification accuracy, the Kappa coefficient, the user accuracy, and the producer accuracy.

[0058] The beneficial effects of the present invention are:

[0059] This application first obtains the original remote sensing image, classifies the original remote sensing image, then creates buffer areas for the divided original remote sensing images respectively, divides the buffer areas into buffer area groups according to the categories of the original remote sensing images, performs spatio-temporal fusion on different buffer area groups through the spatio-temporal fusion algorithm, splices the fused images to generate the final image, and finally performs land use classification using the supervised classification method based on the final image and historical images. In this application, based on the spatio-temporal fusion technology, by integrating the temporal dynamic information of low-resolution images and the spatial details of high-resolution images, an effective solution is provided for improving the accuracy of land use classification. BRIEF DESCRIPTION OF THE DRAWINGS

[0060] In order to more clearly illustrate the technical solutions in the embodiments of the present invention or the prior art, the following will briefly introduce the drawings required for the description of the embodiments or the prior art. Obviously, the drawings in the following description are only some embodiments of the present invention. For those of ordinary skill in the art, other drawings can be obtained based on these drawings without creative efforts.

[0061] Figure 1 is a flowchart provided by an embodiment of a land use classification method based on spatio-temporal fusion technology of the present invention. DETAILED DESCRIPTION OF THE EMBODIMENTS

[0062] To make the objectives, technical solutions, and advantages of the present invention clearer, the technical solutions of the present invention will be described in detail below. Obviously, the described embodiments are only a part of the embodiments of the present invention, rather than all the embodiments. All other implementation manners obtained by those of ordinary skill in the art based on the embodiments of the present invention without creative efforts fall within the scope of protection of the present invention.

[0063] Figure 1 is a flowchart provided by an embodiment of a land use classification method based on spatio-temporal fusion technology of the present invention. Please refer to Figure 1 , this embodiment may include the following steps:

[0064] S1. Obtain the original remote sensing image and classify the original remote sensing image;

[0065] S2. Create buffer areas for the divided original remote sensing images respectively, and divide the buffer areas into buffer area groups according to the categories of the original remote sensing images;

[0066] S3. Perform spatio-temporal fusion on different buffer area groups through the spatio-temporal fusion algorithm;

[0067] S4. Splice the fused images to generate the final image;

[0068] S5. According to the final image and historical images, use the supervised classification method for land use classification.

[0069] It can be understood that in this embodiment, the original remote sensing image is first obtained and classified. Then, buffer areas are created for the divided original remote sensing images respectively, and the buffer areas are divided into buffer area groups according to the categories of the original remote sensing images. The spatio-temporal fusion algorithm is used to perform spatio-temporal fusion on different buffer area groups, the fused images are mosaicked to generate the final image, and finally, according to the final image and historical images, the supervised classification method is used for land use classification. In this embodiment, based on the spatio-temporal fusion technology, by integrating the time dynamic information of the low-resolution image and the spatial details of the high-resolution image, an effective solution is provided for improving the accuracy of land use classification.

[0070] Preferably, step S1 includes:

[0071] Obtain the original remote sensing image;

[0072] Use the supervised classification method to compare the feature vector of each pixel point P of the original remote sensing image with the training sample set to obtain the class label;

[0073] Classify the original remote sensing image according to the class label; wherein, the classification categories include: building area, vegetation area, and other areas;

[0074] Extract the building area and the vegetation area as vector files.

[0075] It can be understood that the original low-resolution remote sensing image is preliminarily classified, and the classification categories include cultivated land, forest land, grassland, building, water body, etc.

[0076] Formula:

[0077] For the classification of each pixel point P, use the supervised classification method, such as support vector machine (SVM), random forest, etc., calculate the feature vector F(p) of this pixel, and compare it with the training sample set to obtain the class label C (p) :

[0078]

[0079] where C is the set of all categories, μ c is the mean feature vector of category c, ‖·‖ 2 is the Euclidean distance.

[0080] Extract the building area and the vegetation area as vector files (such as SHP files).

[0081] Preferably, step S2 includes:

[0082] Dynamically adjust the buffer width of the building area according to the building density of the building area;

[0083] Dynamically adjust the buffer width of the vegetation area according to the vegetation coverage and time dynamics of the vegetation area;

[0084] According to the buffer width of the building area and the buffer width of the vegetation area, use the buffer generation formula to generate the buffer of the building area and the buffer of the vegetation area respectively;

[0085] For the overlapping part of the building area buffer and the vegetation area buffer, retain the information of the high-priority area according to the preset priority, cut off the low-priority area, and introduce a fuzzy allocation mechanism at the junction of the building area buffer and the vegetation area buffer;

[0086] According to the cropped building area buffer and vegetation area buffer, use ArcGIS tools to create corresponding building buffer groups, corresponding vegetation buffer groups and other area boundary groups, and convert them into SHP files.

[0087] It can be understood that the creation of the buffer is based on the characteristics of the ground object type, and the buffer width is dynamically adjusted to adapt to the characteristic requirements of different ground objects. The purpose of the buffer is to expand the scope of the target area to more clearly extract the boundary and optimize the subsequent spatio-temporal fusion processing.

[0088] First, dynamically adjust the buffer width. The dynamic adjustment logic of the buffer width (ddd) is based on the ground object type and regional characteristics, and the rules are as follows:

[0089] Building area buffer:

[0090] Width range: 50m to 100m.

[0091] Basis for dynamic adjustment:

[0092] Building density: Areas with dense buildings require smaller buffers to avoid excessive overlap between buffers; in sparse areas, larger buffers are used to capture more surrounding information.

[0093]

[0094] Among them, ρ is the number of buildings per unit area, ρ max is the maximum value of the building density, k is the adjustment coefficient. β is the boundary complexity: complex boundaries require larger buffers to capture complete boundary information. The boundary complexity is calculated by the gradient change of the building edge

[0095]

[0096] Among them, is the gradient of the building image.

[0097]

[0098] Among them, W is the size of the sliding window (such as 100m × 100m), and N is the number of building pixels within the sliding window

[0099] Vegetation area buffer:

[0100] Width range: 100m to 200m.

[0101] Basis for dynamic adjustment:

[0102] Basis for dynamic adjustment:

[0103] Vegetation coverage (C v ): Larger buffers are used for high-vegetation-coverage areas (such as forest land) to capture the boundary details of the vegetation; smaller buffers are used for low-coverage areas (such as sparse grassland).

[0104] Temporal dynamics (Δ t ): Larger buffers are needed for areas with rapid vegetation growth changes to capture their temporal dynamic characteristics:

[0105]

[0106] Among them, C v is the vegetation coverage ratio per unit area, and C v,max is the maximum vegetation coverage.

[0107] Δ t = NDVI t2 - NDVI t1 , representing the temporal change of the vegetation index of the low-resolution image.

[0108] Then, buffers are generated. The generation formulas for buffers of each type of ground object (buildings, vegetation) are as follows:

[0109] B R = {p|distance(p,R) ≤ d}

[0110] R is the target area (such as a building or vegetation area); d is the width of the buffer after dynamic adjustment; distance(p,R) is the distance from pixel point p to the boundary R of the target area.

[0111] Then image cropping is performed. In the actual processing process, different buffers may overlap (such as the overlap between building buffers and vegetation buffers). To avoid the interference of redundant information on classification and fusion, priority rules for area cropping need to be set.

[0112] The priority is determined according to the significance of the ground features and the requirements of the classification target, and the rules are as follows:

[0113] Priority order: building buffer zone > vegetation buffer zone > other areas.

[0114] Building buffer zone: has the highest priority because of its complex boundary and higher classification requirements; Vegetation buffer zone: the second highest priority, mainly used to capture vegetation dynamics and boundaries; Other areas: the lowest priority, only used for the analysis of background or remaining ground features.

[0115] When buffer zones overlap, the information of the high-priority area is retained according to the priority, and the low-priority area is cropped off.

[0116] The assignment of pixel values in the overlapping part is based on the regional characteristics:

[0117] If a pixel belongs to both the building and vegetation buffer zones at the same time, it is assigned to the building buffer zone.

[0118] For the overlapping area, the assignment formula is:

[0119]

[0120] C(p): the final assigned category of pixel p; R: the set of buffer zones (such as building buffer zone and vegetation buffer zone); P(p|r): the probability that the pixel belongs to region r, determined by the buffer zone priority weight.

[0121] At the junction of the building and vegetation buffer zones (such as rural areas), a fuzzy assignment mechanism can be introduced:

[0122] C(p) = w building ·P building (p) + w vegetation ·P vegetation (p)

[0123] w building ,w vegetation are the assignment weights for the building and vegetation respectively; P building (p) and P vegetation (p) represent the probability values that pixel p belongs to the building and vegetation categories respectively.

[0124] In order to better distinguish the boundaries between the building and vegetation areas and avoid errors in the classification of low-resolution images, it is necessary to expand the building area and vegetation area outwards (the expansion distance is adjusted according to the actual situation) to form buffer zones.

[0125] Finally, create buffers. Using tools such as ArcGIS, create buffers for the building area and the lake area respectively, divide the image area into three categories: building buffer group, vegetation buffer group, and other area boundary group (the remaining part after removing the building buffer and the vegetation buffer), and convert them into SHP files.

[0126] B A ={p|distance(p,A)≤D}

[0127] D is the buffer distance and can be adjusted according to the actual situation.

[0128] Preferably, step S3 includes:

[0129] Use the STDF-A algorithm to perform spatio-temporal fusion on the building buffer group;

[0130] Use the FIT-FIC algorithm to perform spatio-temporal fusion on the vegetation buffer group;

[0131] Use the FSDAF algorithm to perform spatio-temporal fusion on the other area boundary group through a multi-scale feature extraction method.

[0132] It should be noted that FSDAF is a general spatio-temporal fusion algorithm. Based on the multi-scale analysis of the feature space, it can perform fine feature extraction and fusion on images. This algorithm can adapt to various land cover types (such as water bodies, cultivated land, grasslands, etc.) and performs excellently in the overall fusion effect. Its core is to extract image details through multi-scale filtering and then perform adaptive fusion according to the feature weights to ensure the balance of the image in different feature spaces.

[0133] The FIT-FIC algorithm is particularly suitable for spatio-temporal fusion in vegetation-covered areas. The time dynamic characteristics of the vegetation area (such as seasonal changes, fluctuations in vegetation coverage) are significant. This algorithm optimizes the weight allocation of time series information, highlights the spatio-temporal characteristics of vegetation coverage, and at the same time improves the spatial detail performance. Therefore, the fusion effect of FIT-FIC in vegetation areas (such as cultivated land, forest land, grassland) is far better than other algorithms.

[0134] The STDF-A algorithm is designed for the spatio-temporal fusion problem of the building area, focusing on the clear extraction of building boundaries and the restoration of internal details of buildings. The dynamic changes in the building area are usually low, but the boundaries are complex and require high-resolution performance. The STDF-A algorithm combines the edge features of high-resolution images with the time information of low-resolution images and optimizes the fusion effect of the building area through dynamic weight adjustment.

[0135] It can be understood that the low-resolution images and high-resolution images cropped from the regional groups are combined respectively, and appropriate spatio-temporal fusion algorithms are used for each group.

[0136] For the building buffer group:

[0137] Use the STDF-A algorithm (Spatial-Temporal Fusion Algorithm for Building Areas) to fuse low-resolution images and high-resolution images.

[0138] Core formula:

[0139]

[0140] The fused high-resolution image.

[0141] ω L , ω H : Weights of the low-resolution and high-resolution images, satisfying ω L +ω H = 1.

[0142] The feature enhancement function is used to combine the mapping result G(I L ) of the low-resolution image and the details of the high-resolution image.

[0143] For the vegetation buffer group:

[0144] Use the FIT-FIC algorithm to optimize the spatio-temporal characteristics of the vegetation area.

[0145] Core formula:

[0146]

[0147] FIT-FIC optimizes the fusion result through local weighting, and the weights ω L , ω H are dynamically adjusted.

[0148] For other regional boundary groups:

[0149] Use the FSDAF algorithm to fuse water bodies and other regions through a multi-scale feature extraction method.

[0150] Core formula:

[0151]

[0152] Where, F K (I L ) is the k-level filtered feature of the low-resolution image, w k is the weight coefficient, and K is the number of filtering layers.

[0153] Preferably, step S4 includes:

[0154] Determine the reference image and the adjustment image in the fused image;

[0155] Perform gray-scale adjustment and stitching on the adjusted image according to the reference image;

[0156] Perform color normalization on the adjusted image after gray-scale adjustment according to the reference image;

[0157] Perform boundary smoothing on the stitched and fused image to generate the final image.

[0158] It can be understood that color balance and boundary smoothing operations need to be performed before stitching the fused images.

[0159] Color balance: Unify the color characteristics (such as brightness, contrast, hue) of images in different regions to avoid visual inconsistencies caused by fusion algorithms or differences in original images.

[0160] Boundary smoothing: Reduce abrupt transitions in the boundary regions of stitched images in different partitions to achieve seamless connection.

[0161] Preferably, performing gray-scale adjustment and stitching on the adjusted image according to the reference image includes:

[0162] Calculate the gray-scale histograms of the reference image and the adjusted image respectively to determine the gray-scale histogram of the reference image and the gray-scale histogram of the adjusted image;

[0163] Find the cumulative distribution functions of the gray-scale histogram of the reference image and the gray-scale histogram of the adjusted image respectively;

[0164] Perform gray-scale adjustment on the adjusted image according to the cumulative distribution function of the gray-scale histogram of the reference image and the cumulative distribution function of the gray-scale histogram of the adjusted image, and stitch the adjusted adjusted image and the reference image to complete the gray-scale adjustment of the adjusted image.

[0165] It can be understood that color balance aims to adjust the color characteristics of the image to make the images in the stitching area look consistent. Histogram matching is a commonly used color correction method that adjusts the histogram of the image in the area to be stitched to the histogram of the reference image. First, select the area with the largest area or the area with relatively uniform overall brightness in the fused image as the reference, and denote the reference image as I ref , and the image to be adjusted as I target .

[0166] Then calculate the gray-scale histograms H ref and H target of I ref and I target respectively.

[0167] Find the cumulative distribution function (CDF) of the histogram:

[0168]

[0169] where g is the pixel gray value.

[0170] Finally, for each gray value g, find the target gray value g that makes C target (g) closest to C ref (g): new :

[0171] g new = arg min|C target (g) - C ref (g’)|

[0172] Apply g new to I target .

[0173] Concatenate the adjusted I target with I ref .

[0174] Preferably, perform color normalization on the adjusted image after gray value adjustment according to the reference image, including:

[0175] Calculate the mean and standard deviation of the reference image and the adjusted image respectively;

[0176] According to the mean and standard deviation of the reference image, adjust the mean and standard deviation of the adjusted image through the adjustment formula to complete the color normalization of the adjusted image;

[0177] Adjust the brightness and contrast of the concatenated area of the concatenated image.

[0178] It can be understood that color normalization makes the color distribution of the concatenated area consistent by adjusting the mean and standard deviation of the image.

[0179] Calculate the mean μ ref and μ target , μ ref , μ target and the standard deviation σ ref , σ target .

[0180] Complete color normalization according to the adjustment formula. The adjustment formula is:

[0181]

[0182] If the brightness and contrast of the image concatenated area differ greatly, the brightness and contrast can be adjusted separately:

[0183] I target,new (x, y) = α · I target (x, y) + β

[0184] α: Contrast adjustment factor (typical value is 0.8 - 1.2); β: Brightness adjustment factor (typical value is -50 to 50).

[0185] Preferably, perform boundary smoothing on the spliced and fused image to generate the final image, including:

[0186] Define the junction area for the spliced and fused image;

[0187] Perform weight assignment on the junction area;

[0188] Perform weighted summation on the pixel values of the junction area according to the assigned weights to complete the weighted average fusion of the fused image;

[0189] Construct a Gaussian pyramid for the weighted average fused image area, and decompose the image area into sub-images with multiple resolutions;

[0190] Generate a Laplacian pyramid in the Gaussian pyramid to extract the detail information of sub-images with different resolutions;

[0191] Generate a weight pyramid corresponding to the Laplacian pyramid based on the weight calculation of the junction area;

[0192] Perform weighted fusion on the Gaussian pyramid, Laplacian pyramid, and weight pyramid at each resolution layer;

[0193] Reconstruct the fused pyramid from low resolution to high resolution in sequence to generate the final image.

[0194] It can be understood that boundary smoothing aims to eliminate the abrupt transition at the junction of the splicing area, and common methods include weighted average fusion and multi-resolution pyramid fusion.

[0195] Weighted average fusion includes:

[0196] Define the junction area as:

[0197] Extract the junction area B of the spliced image, usually set as a buffer zone with a width of d pixels on both sides of the boundary.

[0198] The weight assignment is:

[0199] Within the junction area, assign weights according to the distance:

[0200] w right (x,y) = 1 - w left (x,y)

[0201] where distance(x,y) represents the distance of the pixel to the boundary, and d is the width of the junction area.

[0202] The weighted average calculation is:

[0203] Perform weighted summation on the pixel values of the boundary region:

[0204] I smooth (x,y) = w left (x,y)·I left (x,y) + w right (x,y)·I right (x,y)

[0205] Multi - resolution pyramid fusion includes:

[0206] Multi - resolution pyramid fusion can handle complex boundaries and is especially suitable for cases where there are large differences in image illumination or texture.

[0207] Generate the Gaussian pyramid as:

[0208] For the spliced image regions I left and I right Construct a Gaussian pyramid and decompose it into sub - images of multiple resolutions.

[0209] Generate the Laplacian pyramid as:

[0210] Generate the Laplacian pyramid from the Gaussian pyramid and extract the detail information of different resolutions.

[0211] Construct the weight pyramid as:

[0212] Generate a weight pyramid corresponding to the Laplacian pyramid based on the weight calculation of the boundary region.

[0213] Fuse the pyramid layers as:

[0214] Perform weighted fusion on each resolution layer:

[0215]

[0216] Reconstruct the fused image:

[0217] Reconstruct the fused pyramid from low resolution to high resolution in sequence to generate the final image.

[0218] Preferably, before step S5, it further includes:

[0219] Obtain a reference image;

[0220] Calculate the peak signal - to - noise ratio, structural similarity, and spectral consistency between the final image and the reference image respectively;

[0221] Evaluate the quality of the final image according to the peak signal - to - noise ratio, structural similarity, and spectral consistency.

[0222] It is understandable that after the generation of the spliced high-resolution image, in order to ensure that the quality of the image meets the actual application requirements, it is necessary to evaluate the accuracy of the fusion effect of the image, and the spatio-temporal fusion quality evaluation is mainly quantified through the following several commonly used indicators.

[0223] Peak Signal-to-Noise Ratio (PSNR)

[0224] PSNR is used to measure the similarity between the fused image and the reference image. The higher the value, the better the image quality. The formula is as follows:

[0225]

[0226] MAX: The maximum possible value of the image pixel value (e.g., 255 represents an 8-bit image).

[0227] MSE: Mean Squared Error, calculation formula:

[0228]

[0229] Where, I fusefd (i) and I ref (i) are the pixel values of the fused image and the reference image respectively, and N is the total number of pixels in the image.

[0230] Structural Similarity (SSIM)

[0231] SSIM measures the similarity between the fused image and the reference image in terms of brightness, contrast, and structure. The value range is [0,1][0,1][0,1], and the value closer to 1 indicates a better fusion effect. The formula is as follows:

[0232]

[0233] In the formula: μ x 、μ y are the averages of the low-resolution optical true image and the fused fine-resolution image at time t2, σ x 、σ y are the image variances of the two, C 1 、C 2 are non-zero constants used to ensure the rationality of the index calculation results. The more similar the overall structures of the two images are, the closer the SSIM value is to 1

[0234] Spectral Angle Mapper (SAM)

[0235] SAM measures the spectral consistency of the fused image, mainly used for multi-spectral or hyperspectral images. The formula is as follows:

[0236]

[0237] The smaller the value, the higher the spectral consistency between the fused image and the reference image.

[0238] Preferably, it further includes:

[0239] Obtain the land use classification result and the ground truth data;

[0240] Calculate the overall classification accuracy of the land use classification result through the confusion matrix method;

[0241] Calculate the Kappa coefficient between the land use classification result and the ground truth data;

[0242] Calculate the user accuracy of the land use classification result and the producer accuracy of the ground truth data;

[0243] Evaluate the land use classification result according to the overall classification accuracy, Kappa coefficient, user accuracy and producer accuracy.

[0244] It can be understood that the classification accuracy evaluation adopts the common confusion matrix method and is evaluated based on the land use classification result of the fused image and the ground truth data (or high-resolution reference image).

[0245] The overall classification accuracy (OA) is:

[0246]

[0247] n ii : The number of pixels of the i-th class on the diagonal of the confusion matrix.

[0248] N: The total number of pixels.

[0249] The Kappa coefficient is:

[0250] The Kappa coefficient measures the consistency between the classification result and the true value data, and its value range is [-1, 1]. The closer the value is to 1, the higher the consistency:

[0251]

[0252] P o : The overall accuracy of the actual classification, which is the same as OA.

[0253] P e : The random agreement accuracy:

[0254]

[0255] Among them, n i. and n .i are the total number of actual values and the total number of predicted values of the i-th class respectively.

[0256] The user accuracy (UA) and the producer accuracy (PA) are as follows:

[0257]

[0258] UA i : The proportion of pixels classified as class i that are actually class i.

[0259] PA i : The proportion of pixels that are actually class i and are classified as class i.

[0260] It can be understood that the same or similar parts in the above embodiments can be referred to each other, and the content not detailed in some embodiments can be seen in the same or similar content in other embodiments.

[0261] It should be noted that in the description of the present invention, the terms "first", "second", etc. are only used for descriptive purposes and cannot be understood as indicating or implying relative importance. In addition, in the description of the present invention, unless otherwise stated, the meaning of "a plurality of" refers to at least two.

[0262] Any process or method description in the flowchart or described in other ways herein can be understood as representing a module, segment, or part of code including one or more executable instructions for implementing a specific logical function or process, and the scope of the preferred embodiments of the present invention includes additional implementations, where the functions can be executed in a substantially simultaneous manner or in the reverse order according to the functions involved, rather than in the order shown or discussed, which should be understood by those skilled in the technical field of the embodiments of the present invention.

[0263] It should be understood that each part of the present invention can be implemented by hardware, software, firmware, or a combination thereof. In the above embodiments, multiple steps or methods can be implemented by software or firmware stored in a memory and executed by a suitable instruction execution system. For example, if implemented by hardware, as in another embodiment, any one or a combination of the following well-known technologies in the art can be used: discrete logic circuits with logic gate circuits for implementing logical functions on data signals, application-specific integrated circuits with appropriate combinational logic gate circuits, programmable gate arrays (PGAs), field-programmable gate arrays (FPGAs), etc.

[0264] Those of ordinary skill in the art of this technology can understand that all or part of the steps carried by the methods in the above embodiments can be completed by instructing relevant hardware through a program, and the program can be stored in a computer-readable storage medium. When the program is executed, it includes one or a combination of the steps of the method embodiments.

[0265] In addition, in each embodiment of the present invention, each functional unit can be integrated into a processing module, or each unit can exist physically alone, or two or more units can be integrated into one module. The above-mentioned integrated module can be implemented in the form of hardware or in the form of a software functional module. When the integrated module is implemented in the form of a software functional module and sold or used as an independent product, it can also be stored in a computer-readable storage medium.

[0266] The above-mentioned storage medium can be a read-only memory, a magnetic disk, an optical disk, or the like.

[0267] In the description of this specification, the description with reference to the terms "one embodiment", "some embodiments", "example", "specific example", or "some examples", etc. means that the specific features, structures, materials, or characteristics described in connection with the embodiment or example are included in at least one embodiment or example of the present invention. In this specification, the schematic representation of the above terms does not necessarily refer to the same embodiment or example. Moreover, the specific features, structures, materials, or characteristics described can be combined in a suitable manner in any one or more embodiments or examples.

[0268] Although the embodiments of the present invention have been shown and described above, it can be understood that the above embodiments are exemplary and should not be construed as limiting the present invention. Those of ordinary skill in the art can make changes, modifications, substitutions, and variations to the above embodiments within the scope of the present invention.

Claims

1. A land use classification method based on spatiotemporal fusion technology, characterized in that: include: Acquiring original remote sensing images, and classifying the original remote sensing images; respectively creating buffer areas for the divided original remote sensing images, and dividing the buffer areas into buffer area groups according to the categories of the original remote sensing images; Performing spatiotemporal fusion on different buffer area groups by using a spatiotemporal fusion algorithm; Stitching the fused images to generate the final image; Based on the final image and historical images, a supervised classification method is used to classify land use.

2. The method according to claim 1, characterized in that: The obtaining of original remote sensing images and classifying the original remote sensing images includes: Obtain original remote sensing images; The supervised classification method is used to compare the feature vector of each pixel point P of the original remote sensing image with the training sample set to obtain a category label; Classifying the original remote sensing image according to the category label; wherein the classification categories include: building area, vegetation area and other areas; The building area and the vegetation area are extracted as vector files.

3. The method according to claim 2, characterized in that The step of creating buffer areas for the divided original remote sensing images respectively and dividing the buffer areas into buffer area groups according to the categories of the original remote sensing images comprises: Dynamically adjust the buffer width of the building area according to the building density of the building area; Dynamically adjusting the buffer width of the vegetation area according to the vegetation coverage and time dynamics of the vegetation area; According to the buffer width of the building area and the buffer width of the vegetation area, a buffer zone of the building area and a buffer zone of the vegetation area are generated respectively by a buffer zone generation formula; For the overlapping part of the building area buffer and the vegetation area buffer, the information of the high priority area is retained according to the preset priority, the low priority area is cut off, and a fuzzy allocation mechanism is introduced for the junction of the building area buffer and the vegetation area buffer; According to the clipped building area buffer and the vegetation area buffer, ArcGIS tools are used to create corresponding building buffer groups, corresponding vegetation buffer groups and other area boundary groups, and convert them into SHP files.

4. The method according to claim 3, characterized in that: The performing spatiotemporal fusion on different buffer area groups by using a spatiotemporal fusion algorithm includes: Using STDF-A algorithm to perform spatiotemporal fusion of the building buffer group; Using the FIT-FIC algorithm to perform spatiotemporal fusion of the vegetation buffer group; The FSDAF algorithm is used to perform spatiotemporal fusion on the other region boundary groups through a multi-scale feature extraction method.

5. The method according to claim 4, characterized in that The step of stitching the fused images to generate a final image includes: determining reference images and adjusting images in fused images; Performing grayscale adjustment and splicing on the adjustment image according to the reference image; Performing color normalization on the adjusted image after grayscale adjustment according to the reference image; The stitched and fused images are subjected to boundary smoothing processing to generate a final image.

6. The method according to claim 5, characterized in that The grayscale adjustment and splicing of the adjustment image according to the reference image includes: Calculating grayscale histograms for the reference image and the adjusted image respectively, and determining a reference image grayscale histogram and an adjusted image grayscale histogram; Calculating cumulative distribution functions of the reference image grayscale histogram and the adjusted image grayscale histogram respectively; Grayscale adjustment of the adjustment image is performed according to the integral distribution function of the grayscale histogram of the reference image and the cumulative distribution function of the grayscale histogram of the adjustment image, and the adjusted adjustment image and the reference image are spliced ​​to complete the grayscale adjustment of the adjustment image.

7. The method according to claim 6, characterized in that The step of performing color normalization on the grayscale-adjusted adjusted image according to the reference image comprises: Calculating the mean and standard deviation of the reference image and the adjusted image respectively; According to the mean and standard deviation of the reference image, the mean and standard deviation of the adjustment image are adjusted by an adjustment formula to complete color normalization of the adjustment image; Adjust the brightness and contrast of the stitched area of ​​the stitched image.

8. The method according to claim 7, characterized in that The step of performing boundary smoothing on the stitched and fused images to generate a final image includes: defining a boundary area for the stitched and fused images; Allocating weights to the boundary area; Performing weighted summation on the pixel values ​​of the boundary area according to the assigned weights to complete weighted average fusion of the fused image; Constructing a Gaussian pyramid for the image area after weighted average fusion, and decomposing the image area into sub-images with multiple layers of resolution; Generating a Laplacian pyramid in the Gaussian pyramid to extract detail information of sub-images with different resolutions; Generate a weight pyramid corresponding to the Laplacian pyramid based on the weight calculation of the boundary area; Performing weighted fusion on the Gaussian pyramid, the Laplacian pyramid, and the weight pyramid at each resolution layer; The fused pyramids are reconstructed from low resolution to high resolution to generate the final image.

9. The method according to claim 8, characterized in that Before performing land use classification using a supervised classification method based on the final image and the historical image, the method further includes: Obtain reference images; respectively calculating the peak signal-to-noise ratio, structural similarity and spectral consistency of the final image and the reference image; The quality of the final image is evaluated according to the peak signal-to-noise ratio, the structural similarity and the spectral consistency.

10. The method according to claim 9, characterized in that Also includes: Obtain land use classification results and ground truth data; Calculating the total classification accuracy of the land use classification result by using a confusion matrix method; Calculating the Kappa coefficient between the land use classification result and the ground truth data; calculating a user precision of the land use classification result and a producer precision of the ground truth data; The land use classification result is evaluated according to the overall classification accuracy, the Kappa coefficient, the user accuracy and the producer accuracy.

Images