A land use classification method and system based on remote sensing technology

Abstract

The present invention relates to the field of land classification and discloses a land use classification method and system based on remote sensing technology, which is used to solve the problem that the illumination angle changes and does not match the boundary of the splicing area when performing land use classification. The method and system include collecting baseline data of sub-detection areas, obtaining real-time multiple remote sensing data for each sub-detection area, evaluating to obtain a comprehensive illumination change index, and judging the illumination angle change. If it is judged that the illumination angle has changed, the image illumination normalization and shadow removal are performed to obtain the sub-detection area image, and the image of each sub-detection area is stitched to obtain an initial stitched image. The stitching error index is evaluated and the stitching consistency is judged according to the stitching error index. If it is judged that the initial stitched image is inconsistent, the stitching optimization of the initial stitched image is performed to obtain an actual stitched image, thereby effectively improving the accuracy of remote sensing image acquisition and improving the rationality of land use classification.

Description

Technical Field

[0001] The present invention relates to the field of land classification, and more particularly to a land use classification method and system based on remote sensing technology. Background Art

[0002] In the field of modern land use classification, remote sensing technology has become an essential tool for analyzing land cover types and changes. Remote sensing imagery enables researchers and decision makers to efficiently perform tasks such as land classification. Land use classification using remote sensing imagery often involves analyzing multi-temporal data to identify dynamic changes in land types and provide accurate classification results.

[0003] Existing remote sensing land use classification methods typically rely on the acquisition of multi-source data, including optical satellite imagery, radar imagery, and elevation imagery, combined with machine learning or deep learning models for classification. The general process includes data preprocessing, feature extraction, classification algorithm training and prediction, and subsequent classification accuracy evaluation. In practical applications, remote sensing images are typically acquired at fixed intervals. After splicing, alignment, and classification, these images can be used for land use type analysis and long-term monitoring.

[0004] However, in the process of implementing the technical solutions of the invention in the embodiments of the present application, the present application found that the above technology has at least the following technical problems:

[0005] In practical applications, even if images are collected at fixed intervals, remote sensing images taken at different times may still display inconsistent representations of the same area in different images due to factors such as lighting conditions and changes in observation angle. For example, due to changes in the solar altitude, the image brightness and shadow length of the same feature at different times may vary, affecting the stability of the classification algorithm. Furthermore, slight shifts in the satellite orbit or differences in imaging angles can lead to mismatched boundaries in the stitched areas. Such errors not only affect the consistency of land use classification results but can also lead to misjudgment of land cover types, thus affecting subsequent analysis and decision-making. Summary of the Invention

[0006] In order to overcome the above-mentioned defects of the prior art, the present invention provides a land use classification method and system based on remote sensing technology to solve the problems existing in the above-mentioned background technology.

[0007] To achieve the above object, the present invention provides the following technical solutions:

[0008] A land use classification method based on remote sensing technology includes the following steps: Step 1: Divide the target area into multiple sub-detection areas, select a benchmark time point, collect multiple remote sensing data for each sub-detection area at the benchmark time point, and use the collected multiple remote sensing data as the benchmark data of the sub-detection area, the remote sensing data includes optical satellite images, radar images and elevation images, and the benchmark data includes benchmark optical satellite images, benchmark radar images and benchmark elevation images; Step 2: Set a fixed collection time period, for each sub-detection area, obtain real-time multiple remote sensing data at the current collection time point, the real-time multiple remote sensing data includes real-time optical satellite images and real-time radar images; Step 3: Evaluate the comprehensive light change index of each sub-detection area based on the real-time multiple remote sensing data and the benchmark data, and use the comprehensive light change index to obtain the real-time multiple remote sensing data. The illumination change index is used to judge the change of illumination angle; Step 4: If it is judged that the illumination angle has changed, the real-time optical satellite image is subjected to image illumination normalization and shadow removal to obtain a corrected optical satellite image, and the sub-detection area image is obtained based on the corrected optical satellite image, combined with the radar image and the elevation image; Step 5: Based on the same acquisition time point, the images of each sub-detection area are stitched to obtain an initial stitched image, and the stitching error index is obtained based on the evaluation of the initial stitched image. The stitching consistency of the initial stitched image is judged based on the stitching error index; Step 6: If it is judged that the initial stitched image is inconsistent, the initial stitched image is stitched and optimized to obtain the actual stitched image; Step 7: The actual stitched image of each acquisition time point is stored in the database, and land use classification is performed in combination with the deep learning method.

[0009] Preferably, the steps for obtaining the comprehensive illumination change index are: obtaining main band data of the real-time optical satellite image and the reference optical satellite image, the main band data including the red, short-wave infrared and near-infrared bands of each pixel, and evaluating the optical image change coefficient based on the main band data of the real-time optical satellite image and the reference optical satellite image; obtaining the polarization component of each pixel of the real-time radar image and the reference radar image, the polarization component including the vertical-vertical polarization component and the vertical-horizontal polarization component, and calculating the radar change coefficient based on the polarization component of each pixel; obtaining the reference elevation image, and evaluating the elevation influence coefficient based on the reference elevation image; normalizing the optical image change coefficient, the radar change coefficient and the elevation influence coefficient, and evaluating the normalized optical image change coefficient, the radar change coefficient and the elevation influence coefficient to obtain the comprehensive illumination change index. The specific acquisition steps are: Where LC is the comprehensive illumination change index, IC is the optical image change coefficient, RCC is the radar change coefficient, EIC is the elevation influence coefficient, and a1, a2, and a3 are the weight coefficients of the optical image change coefficient, the radar change coefficient, and the elevation influence coefficient.

[0010] Preferably, the optical image variation coefficient acquisition step is as follows: for each pixel point, calculating the spectral change between the current time point and the reference time point, the spectral change including the red light difference, the short-wave infrared difference and the near-infrared difference; taking the spectral change of the pixel point as the feature vector, taking the spectral change of each pixel point as the data set, and the spectral change of each pixel point in the data set as the data point; using the silhouette coefficient method to obtain the optimal number of clusters, using the K-means clustering method to cluster the data set to obtain the final cluster cluster; calculating the spectral change mean of each cluster cluster based on the final cluster cluster, and calculating the optical image variation coefficient based on the spectral change mean of each cluster cluster. The specific acquisition steps are as follows: Where IC is the optical image variation coefficient, K is the number of clusters, and C k It is expressed as the mean spectral change of the kth cluster.

[0011] Preferably, the steps of clustering the data set using the K-means clustering method to obtain the final cluster cluster are as follows: Step 3.1: randomly select K data points in the data set as initial cluster centers, and for each data point, calculate its Euclidean distance to each initial cluster center. For each data point, traverse the K initial cluster centers and assign it to the cluster corresponding to the nearest initial cluster center; Step 3.2: traverse all data points to obtain initial cluster clusters, and for each initial cluster cluster, calculate the mean of the data points within it to obtain a new cluster center; Step 3.3: repeat steps 3.1 and 3.2 until the cluster center no longer changes, and obtain the final cluster cluster.

[0012] Preferably, the radar variation coefficient acquisition step comprises: acquiring a vertical-vertical polarization component and a vertical-horizontal polarization component of each pixel point of the real-time radar image and the reference radar image;

[0013] For each pixel, the logarithmic ratio change rate is calculated based on the vertical-vertical polarization component and the vertical-horizontal polarization component. The radar variation coefficient is calculated based on the logarithmic ratio change rate of the vertical-vertical polarization component and the vertical-horizontal polarization component of each pixel. The specific acquisition steps are as follows: Where RCC is the radar variation coefficient, N is the number of pixels, Δσ 0,i Expressed as the logarithmic ratio change rate of the vertical-vertical polarization component of the i-th pixel, Δσ 1,i It is expressed as the logarithmic ratio change rate of the vertical to horizontal polarization components of the i-th pixel.

[0014] Preferably, the elevation influence coefficient acquisition step is as follows: for each pixel point in the reference elevation image, an elevation value is acquired, and the slope and slope direction of each pixel point are calculated based on the elevation value of each pixel point; a reference solar altitude angle and a reference solar azimuth angle are acquired based on the reference optical satellite image, and a real-time solar altitude angle and a real-time solar azimuth angle are acquired based on the real-time optical satellite image; the solar incident angle at the reference time point is calculated based on the reference solar altitude angle and the reference solar azimuth angle, and the real-time solar incident angle is calculated based on the real-time solar altitude angle and the real-time solar azimuth angle; the elevation influence coefficient is calculated based on the solar incident angle. The specific acquisition steps are as follows: Where EIC is the elevation influence coefficient, N is the number of pixels, and θ i (T ref ) is the solar incident angle of the i-th pixel at the reference time point, θ i (T) is the real-time solar incident angle of the i-th pixel.

[0015] Preferably, the step of judging the change of illumination angle based on the comprehensive illumination change index is: comparing the comprehensive illumination change index with the illumination change threshold; if the comprehensive illumination change index is less than the illumination change threshold, it is judged that the illumination angle has not changed; if the comprehensive illumination change index is greater than or equal to the illumination change threshold, it is judged that the illumination angle has changed.

[0016] Preferably, the stitching error index obtaining step comprises: stitching the images of each sub-detection area to obtain an initial stitched image, obtaining the image overlap area in the initial stitched image, obtaining the spectral value of each pixel point of the two sub-detection areas in the image overlap area through the optical satellite images of the sub-detection areas, and calculating the average spectral difference of all image overlap areas in the initial stitched image; obtaining the brightness value of each pixel point of the two sub-detection areas in the image overlap area through the optical satellite images of the sub-detection areas, and calculating the average brightness difference of all image overlap areas in the initial stitched image; and obtaining the average brightness difference of all image overlap areas in the initial stitched image through the Harris The corner detection method extracts key points in the overlapping area of ​​the images of each sub-detection area, matches the key points of the two overlapping sub-detection areas using the FLANN method, calculates the Euclidean distance of the matching points using the Euclidean distance, and averages the Euclidean distances of all matching points in the initial stitched image to obtain the average geometric offset error; the average spectral difference, average brightness difference, and average geometric offset error are normalized, and the stitching error index is calculated based on the normalized average spectral difference, average brightness difference, and average geometric offset error. The specific acquisition steps are: SEI = ln(1+SD)+BD+e GO ; Where SEI is the stitching error index, SD is the average spectral difference, BD is the average brightness difference, and GO is the average geometric offset error.

[0017] Preferably, the step of judging the stitching consistency of the initial stitched image according to the stitching error index is: comparing the stitching error index with the error threshold; if the stitching error index is less than the error threshold, then judging that the stitching of the initial stitched image is consistent; if the stitching error index is greater than or equal to the error threshold, then judging that the stitching of the initial stitched image is inconsistent.

[0018] Preferably, a land use classification system based on remote sensing technology includes: a benchmark data acquisition module, which divides the target area into multiple sub-detection areas, selects a benchmark time point, collects multiple remote sensing data for each sub-detection area at the benchmark time point, and uses the collected multiple remote sensing data as the benchmark data of the sub-detection area, the remote sensing data includes optical satellite images, radar images and elevation images, and transmits the benchmark data to the illumination angle change judgment module; a multiple remote sensing data acquisition module is used to set a fixed acquisition time period, and for each sub-detection area, acquires real-time multiple remote sensing data at the current acquisition time point, the real-time multiple remote sensing data includes real-time optical satellite images and real-time radar images, and transmits the real-time multiple remote sensing data to the illumination angle change judgment module; the illumination angle change judgment module is used to evaluate the comprehensive illumination change index of each sub-detection area based on the real-time multiple remote sensing data and the benchmark data, and perform illumination angle evaluation based on the comprehensive illumination change index. Change judgment; illumination correction module, if it is judged that the illumination angle has changed, the real-time optical satellite image is subjected to image illumination normalization and shadow removal to obtain a corrected optical satellite image, and the sub-detection area image is obtained based on the corrected optical satellite image in combination with the radar image and the elevation image, and the sub-detection area image is transmitted to the stitching consistency judgment module; the stitching consistency judgment module is used to stitch each sub-detection area image based on the same acquisition time point to obtain an initial stitched image, and obtain a stitching error index based on the initial stitched image evaluation, and perform stitching consistency judgment on the initial stitched image based on the stitching error index; the stitching optimization module, if it is judged that the initial stitched image is inconsistent, the initial stitched image is stitched and optimized to obtain an actual stitched image, and the actual stitched image is transmitted to the land use classification module; the land use classification module is used to store the actual stitched images of each acquisition time point in the database, and perform land use classification in combination with the deep learning method.

[0019] Technical effects and advantages of the present invention:

[0020] Collect benchmark data of sub-detection areas, obtain real-time multiple remote sensing data for each sub-detection area, evaluate and obtain a comprehensive illumination change index, and judge the illumination angle change. If it is judged that the illumination angle has changed, perform image illumination normalization and shadow removal to obtain the sub-detection area image, stitch each sub-detection area image to obtain an initial stitched image, evaluate and obtain a stitching error index, and judge the stitching consistency based on the stitching error index. If it is judged that the initial stitched image is inconsistent, perform stitching optimization on the initial stitched image to obtain the actual stitched image, effectively improving the accuracy of remote sensing image acquisition and the rationality of land use classification. BRIEF DESCRIPTION OF THE DRAWINGS

[0021] Figure 1 Flowchart of a land use classification method based on remote sensing technology provided in this application embodiment

[0022] Figure 2 A structural diagram of a land use classification system based on remote sensing technology provided in an embodiment of the present application. DETAILED DESCRIPTION

[0023] The technical solutions of the present invention will be clearly and completely described below in conjunction with the drawings in the present invention. In addition, the forms of the various structures described in the following embodiments are merely examples. The land use classification method and system based on remote sensing technology involved in the present invention are not limited to the various structures described in the following embodiments. All other implementations obtained by ordinary technicians in this field without making creative work are within the scope of protection of the present invention.

[0024] The present invention provides a land use classification method based on remote sensing technology, such as Figure 1 As shown, the following steps are included:

[0025] Step 1: Divide the target area into multiple sub-detection areas, select a benchmark time point, and collect multiple remote sensing data for each sub-detection area at the benchmark time point. The collected multiple remote sensing data are used as the benchmark data of the sub-detection area for subsequent comparative analysis. The remote sensing data includes optical satellite images, radar images, and elevation images, and the benchmark data includes benchmark optical satellite images, benchmark radar images, and benchmark elevation images.

[0026] When the target area is divided into multiple sub-detection areas, a division method can be selected according to actual conditions, including division based on fixed grids, division based on administrative regions, and division based on natural geographical units;

[0027] Fixed grid division is based on fixed longitude and latitude grids, with common grid sizes being 1km×1km, 10km×10km, and 100km×100km; administrative division is based on administrative units (such as provinces, cities, counties, or townships); and natural geographical unit division is based on natural geographical features such as watershed distribution, mountains, plains, and wetlands.

[0028] In remote sensing image analysis, selecting a benchmark time point is crucial because it directly affects the accuracy of subsequent illumination normalization, image stitching, and land use classification. The benchmark time point should be when image quality is optimal, lighting conditions are stable, and cloud interference is minimized.

[0029] In this embodiment, it should be specifically explained that the steps of selecting the reference time point are:

[0030] Identify the key factors that affect image quality, including lighting conditions (moderate solar altitude), meteorological conditions (no or low cloud cover), ground stability (such as vegetation growth and urban expansion), and the availability of remote sensing data (optical satellite imagery, radar imagery, and elevation imagery are all available). Select a time point that meets all of these key factors as the initial reference time point.

[0031] Acquire multi-temporal remote sensing data for the target area over the past two to three years. This data includes optical satellite imagery, radar imagery, elevation imagery, and meteorological data, which are used to analyze changes in illumination, cloud cover, and surface features at different points in time.

[0032] Cloud cover can significantly impact the usability of remote sensing images, so it's important to first select images with minimal cloud cover. Scene classification maps for optical satellite imagery can be used to automatically detect cloud cover and eliminate images with cloud cover exceeding 10%. At the same time, candidate images are manually inspected to ensure they are free of light cloud, haze, or atmospheric scattering to ensure clarity and usability.

[0033] The sun's altitude directly affects the brightness and shadow range of the image, so it is necessary to select time points with moderate sun angles. Sun angle information can be extracted from remote sensing image metadata, and images with sun angles between 40° and 70° can be selected to reduce the effects of long shadows and overexposure. Furthermore, a curve of the sun's altitude variation throughout the year can be plotted to ensure stable lighting conditions at the selected time points.

[0034] The stability of surface features is crucial for illumination normalization and subsequent classification. Therefore, key indicators such as vegetation index, building index, and water index need to be calculated to select the time point with the least surface change. For example, the stability of vegetation growth can be calculated through vegetation index, and the time point with the smallest vegetation index change can be selected to reduce spectral variation errors caused by factors such as crop growth and seasonal changes in forests.

[0035] Check whether there is a corresponding radar image at the initial benchmark time point for shadow detection and lighting correction. In addition, although elevation images are usually fixed data, it is still necessary to ensure that the resolution and quality can support image analysis to avoid inaccurate terrain correction due to elevation image errors;

[0036] Taking all screening factors into consideration, a time point was selected from among all initial benchmark time points to meet the following criteria: no or low cloud cover, moderate solar altitude, stable surface reflectivity, and the availability of optical, radar, and elevation imagery. Finally, image quality was manually reviewed to ensure the benchmark imagery was free of human interference (such as road construction and pollution). This time point served as the benchmark time point for the sub-inspection area and was used for subsequent image stitching, illumination normalization, and land use classification analysis.

[0037] Step 2: Set a fixed acquisition time period. For each sub-detection area, obtain real-time multiple remote sensing data at the current acquisition time point. The real-time multiple remote sensing data includes real-time optical satellite images and real-time radar images.

[0038] Step 3: Based on the real-time multiple remote sensing data and benchmark data, the comprehensive illumination change index of each sub-detection area is obtained, and the illumination angle change is judged based on the comprehensive illumination change index;

[0039] In this embodiment, it should be specifically explained that the steps for obtaining the comprehensive illumination change index are:

[0040] Obtain the main band data of real-time optical satellite images and benchmark optical satellite images. The main band data includes the red, short-wave infrared, and near-infrared bands of each pixel point. The optical image variation coefficient is evaluated based on the main band data of the real-time optical satellite images and the benchmark optical satellite images.

[0041] Obtain the polarization component of each pixel in the real-time radar image and the baseline radar image. The polarization component includes vertical-vertical polarization component and vertical-horizontal polarization component. Calculate the radar variation coefficient based on the polarization component of each pixel. The radar effect is not affected by illumination, but is affected by changes in surface humidity and terrain roughness. Therefore, it can be used to help determine which changes are due to illumination and which are due to changes in ground objects.

[0042] Obtain benchmark elevation images and evaluate the elevation influence coefficient based on the benchmark elevation images;

[0043] Real-time multiple remote sensing data do not include elevation images because elevation images usually do not change over time. Therefore, when calculating the elevation influence coefficient, there is no need to compare the baseline elevation image and the current elevation image. Traditional change detection tasks (such as terrain change monitoring) may calculate the difference between the baseline elevation image and the current elevation image, but in illumination normalization, we only need to care about how the terrain affects the incident angle of sunlight, rather than whether the terrain itself changes. Since the solar altitude angle and solar azimuth angle change over time, areas facing the sun will receive more light, while areas facing away from the sun will darken or form shadows. By calculating the slope and aspect through a fixed elevation image, we can infer the solar incidence angle at different time points and quantify the impact of the terrain on illumination changes.

[0044] The optical image variation coefficient, radar variation coefficient, and elevation influence coefficient are normalized, and the comprehensive illumination variation index is obtained based on the normalized optical image variation coefficient, radar variation coefficient, and elevation influence coefficient. The specific acquisition steps are as follows:

[0045]

[0046] Where LC represents the comprehensive illumination variation index, and IC represents the optical image variation coefficient. A large optical image variation coefficient indicates significant changes in image characteristics such as brightness, color, and shadows, which may be caused by changes in sun angle, weather conditions, or surface reflectivity. This change directly affects the comprehensive illumination variation index, resulting in significant variations in illumination conditions at different points in time, thereby affecting the accuracy of remote sensing image illumination normalization and land use classification. RCC represents the radar variation coefficient. Since radar images are not affected by illumination, their variations are primarily due to changes in surface moisture, roughness, or surface structure. A high radar variation coefficient indicates that surface variation, rather than illumination, is the primary factor, and the comprehensive illumination variation index may be low. A low radar variation coefficient indicates that the surface is relatively stable and illumination variation is the primary factor, and the comprehensive illumination variation index may be high. EIC represents the elevation influence coefficient. Since the sun's altitude and azimuth angles vary at different times, areas with high terrain elevation can cause significant variations in illumination distribution, affecting image characteristics such as brightness and shadow location. When the elevation influence coefficient is high, it means that the changes in slope and aspect have a greater impact on the solar incidence angle, causing significant changes in lighting conditions and an increase in the comprehensive lighting change index. On the contrary, when the terrain is relatively flat, the lighting conditions are relatively stable and the degree of illumination change in the image is small. a1, a2, and a3 are the weight coefficients of the optical image variation coefficient, the radar variation coefficient, and the elevation influence coefficient, and a1+a2+a3=1. a1, a2, and a3 are obtained through the hierarchical analysis method. For example, a1, a2, and a3 can be 0.4, 0.3, and 0.3.

[0047] The Analytic Hierarchy Process (AHP) is a decision analysis method that breaks down complex problems into multiple factors through a hierarchical structure. The weights are then calculated through pairwise comparisons to determine the relative importance of each factor to the final decision. A judgment matrix is ​​constructed to compare the importance of each factor pairwise. Expert or data analysis methods can be used to assign values. Then, through consistency checks, the relative weights of each factor are calculated to ensure that the weights are reasonable and reliable.

[0048] Polarization components are the polarization characteristics of electromagnetic waves as they return to the radar receiver after interacting with the surface. These components primarily include vertical-vertical polarization and vertical-horizontal polarization. Vertical-vertical polarization indicates that both the electromagnetic waves transmitted and received by the radar are vertically polarized and is commonly used for detecting smooth surfaces such as water bodies and urban buildings. Vertical-horizontal polarization indicates that the electromagnetic waves transmitted by the radar are vertically polarized but received horizontally, making it suitable for identifying complex surface features such as vegetation and forests. Analyzing changes in these polarization components can help distinguish between changes in illumination, surface moisture, or changes in the ground objects themselves, improving the classification accuracy of remote sensing images.

[0049] In this embodiment, it should be specifically explained that the steps for obtaining the optical image variation coefficient are:

[0050] For each pixel, calculate the spectral change between the current time point and the reference time point. The spectral change includes the red light difference, short-wave infrared difference, and near-infrared difference.

[0051] The spectral change of the pixel point is used as the feature vector, the spectral change of each pixel point is used as the data set, and the spectral change of each pixel point in the data set is used as the data point;

[0052] The silhouette coefficient method is used to obtain the optimal number of clusters, and the K-means clustering method is used to cluster the data set to obtain the final clusters;

[0053] The spectral variation mean of each cluster is calculated based on the final cluster, and the optical image variation coefficient is calculated based on the spectral variation mean of each cluster. The specific acquisition steps are as follows:

[0054]

[0055] Where IC is the optical image variation coefficient, K is the number of clusters, and C k It is expressed as the mean spectral change of the kth cluster.

[0056] The Silhouette Coefficient method is a technique used to evaluate clustering effectiveness and determine the optimal number of clusters. It determines the rationality of clustering by measuring the closeness of data points within their clusters and their separation from the nearest cluster. This method effectively avoids the degradation of clustering quality caused by too many or too few clusters, making data classification more reasonable and improving the accuracy of subsequent analysis.

[0057] In this embodiment, it should be specifically explained that the steps of clustering the data set using the K-means clustering method to obtain the final clusters are as follows:

[0058] Step 3.1: Randomly select K data points in the data set as the initial cluster centers. For each data point, calculate the Euclidean distance from each initial cluster center. The specific method is to obtain Where d(P,β) represents the Euclidean distance from the data point to the cluster center, where P represents the data point and β represents the initial cluster center. For each data point, traverse the K initial cluster centers and assign it to the cluster corresponding to the nearest initial cluster center.

[0059] Step 3.2: After traversing all data points, we get the initial clusters. For each initial cluster, we calculate the mean of the data points in it and get the new cluster center.

[0060] Step 3.3: Repeat steps 3.1 and 3.2 until the cluster center no longer changes, and the final cluster is obtained.

[0061] In this embodiment, it should be specifically explained that the steps for obtaining the radar variation coefficient are:

[0062] Obtain the vertical-vertical polarization component and vertical-horizontal polarization component of each pixel of the real-time radar image and the reference radar image;

[0063] For each pixel, the logarithmic ratio change rate is calculated based on the vertical-vertical polarization component and the vertical-horizontal polarization component. The specific acquisition steps are as follows:

[0064] Δσ0=|log(σ0 VV (T))-log(σ0 VV (T ref ))|;

[0065] Δσ1=|log(σ0 VH (T))-log(σ0 VH (T ref ))|;

[0066] Where Δσ0 is the rate of change of the logarithmic ratio of the vertical-vertical polarization component, Δσ1 is the rate of change of the logarithmic ratio of the vertical-horizontal polarization component, and σ0 is VV(T) represents the vertical-vertical polarization component of the real-time radar image, σ0 VV (T ref ) is expressed as the vertical-vertical polarization component of the reference radar image, σ0 VH (T) represents the vertical-horizontal polarization component of the real-time radar image, σ0 VH (T ref ) is expressed as the vertical-horizontal polarization component of the reference radar image;

[0067] The radar variation coefficient is calculated based on the logarithmic ratio change rate of the vertical-vertical polarization component and the vertical-horizontal polarization component of each pixel. The specific acquisition steps are as follows:

[0068]

[0069] Where RCC is the radar variation coefficient. If RCC is small, for example, RCC < 0.1, it means that most of the changes come from illumination rather than ground objects. If RCC is large, it means that the ground objects themselves have changed and the influence of illumination is small. N is the number of pixels, and Δσ is the distance between the pixels and the ground objects. 0,i Expressed as the logarithmic ratio change rate of the vertical-vertical polarization component of the i-th pixel, Δσ 1,i It is expressed as the logarithmic ratio change rate of the vertical to horizontal polarization components of the i-th pixel.

[0070] In this embodiment, it should be specifically explained that the steps for obtaining the elevation influence coefficient are:

[0071] For each pixel point in the benchmark elevation image, obtain the elevation value, and calculate the slope and slope direction of each pixel point based on the elevation value of each pixel point;

[0072] Obtain the benchmark solar altitude and azimuth angles based on the benchmark optical satellite imagery, and obtain the real-time solar altitude and azimuth angles based on the real-time optical satellite imagery. The solar altitude is the angle between the sun's rays and the horizon, which determines the intensity of illumination and the length of shadows. The solar azimuth is the angle of the sun relative to true north, which determines the direction of illumination and the direction of shadow projection.

[0073] The solar incident angle at the reference time point is calculated based on the reference solar altitude angle and the reference solar azimuth angle. The real-time solar incident angle is calculated based on the real-time solar altitude angle and the real-time solar azimuth angle. The specific acquisition steps are as follows:

[0074] cos(θ i (T ref ))=cos(Z ref )cos(S i )+sin(Z ref )sin(S i )cos(Aref -A i );

[0075] cos(θ i (T))=cos(Z T )cos(S i )+sin(Z T )sin(S i )cos(A T -A i );

[0076] Where θ i (T ref ) is the solar incident angle of the i-th pixel at the reference time point, θ i (T) is the real-time solar incident angle at the i-th pixel, Z ref With A ref are the reference solar altitude angle and the reference solar azimuth angle, Z T With A T They are the real-time solar altitude angle and the real-time solar azimuth angle. When calculating the solar incidence angle, the slope and slope direction remain unchanged because the terrain remains unchanged.

[0077] The elevation influence coefficient is calculated based on the solar incidence angle. The specific steps are as follows:

[0078]

[0079] Where EIC is the elevation influence coefficient, N is the number of pixels, and the difference in solar incidence angles for all pixels is calculated and averaged to measure the impact of terrain on illumination changes.

[0080] In this embodiment, it should be specifically explained that the steps for calculating the slope and slope direction of each pixel point according to the elevation value of each pixel point are as follows:

[0081] Get the elevation values ​​of the pixel's east-west adjacent pixels and calculate the east-west slope of the pixel. The specific steps are as follows:

[0082]

[0083] Where G x It is expressed as the east-west slope, d is the pixel resolution, H(x+1,y) and H(x-1,y) are the elevation values ​​of the east and west adjacent pixels of the pixel respectively;

[0084] Get the elevation values ​​of the pixel points adjacent to the north and south of the pixel point, and calculate the north-south slope of the pixel point. The specific steps are as follows:

[0085]

[0086] Where G y It is expressed as the north-south slope, d is the pixel resolution, H(x,y+1) and H(x,y+1) are the elevation values ​​of the adjacent pixels to the north and south of the pixel, respectively;

[0087] The slope is calculated based on the east-west slope and the west-east slope. The specific steps are as follows:

[0088]

[0089] Where S is the slope, G x Expressed as the east-west slope, G y It is expressed as north-south slope;

[0090] For each pixel point, the slope direction is calculated based on the east-west slope and the north-south slope. The specific steps for obtaining the slope direction are as follows:

[0091]

[0092] Where A is the slope aspect, which indicates the direction of the slope relative to the north direction and ranges from 0° to 360°. It is used to determine the impact of the sunlight incident angle.

[0093] In this embodiment, it should be specifically explained that the steps for determining the change in illumination angle according to the comprehensive illumination change index are as follows:

[0094] The integrated illumination change index is compared with the illumination change threshold. If the integrated illumination change index is less than the illumination change threshold, the illumination angle is considered unchanged. If the integrated illumination change index is greater than or equal to the illumination change threshold, the illumination angle is considered to have changed. The illumination change threshold is determined using the adaptive thresholding method, an algorithm that dynamically determines the threshold based on data distribution. This method automatically adapts to changes in different scenarios, making the threshold more accurate and stable. When calculating the illumination change threshold, the adaptive thresholding method automatically selects the optimal threshold based on the data distribution of the integrated illumination change index.

[0095] Step 4: If it is determined that the illumination angle has changed, the real-time optical satellite image is subjected to image illumination normalization and shadow removal to obtain a corrected optical satellite image. The corrected optical satellite image is then combined with the radar image and the elevation image to obtain an image of the sub-detection area. It should be noted that image illumination normalization and shadow removal are conventional techniques, and the specific steps are not described in detail in this embodiment.

[0096] Image illumination normalization and shadow removal is a processing method for radiation correction and shadow compensation of optical remote sensing images. It aims to reduce image brightness, contrast and color deviations caused by illumination changes and terrain shadows, so that images have consistent lighting conditions at different times and in different areas, thereby improving the accuracy of image analysis and classification.

[0097] Illumination normalization usually adjusts the image brightness and color to make it close to the reference image through methods such as histogram matching, least squares regression, and terrain illumination correction based on DEM; shadow removal often uses methods such as shadow compensation based on SAR images, spectral mixing analysis, and deep learning restoration to fill in the information of the shadow area and reduce the errors caused by uneven lighting.

[0098] In this embodiment, it should be specifically explained that the steps for obtaining the sub-detection area image based on the corrected optical satellite image, combined with the radar image and the elevation image are as follows:

[0099] The corrected optical satellite images, radar images and elevation images are aligned through geometric correction; the ground feature information is enhanced by combining radar images, and the terrain information is extracted by combining elevation images. The relationship between optical satellite images, radar images and elevation images is learned through neural networks to generate multi-source fused images.

[0100] Step 5: Based on the same acquisition time point, stitch the images of each sub-detection area to obtain an initial stitched image. Evaluate the initial stitched image to obtain a stitching error index, and judge the stitching consistency of the initial stitched image based on the stitching error index.

[0101] In this embodiment, it should be specifically explained that the steps for obtaining the splicing error index are:

[0102] Perform image stitching on each sub-detection area to obtain an initial stitched image, obtain an image overlap area in the initial stitched image, where the image overlap area is the overlapping portion of the images of the two sub-detection areas, obtain the spectral value of each pixel point of the two sub-detection areas in the image overlap area through the optical satellite images of the sub-detection areas, and calculate the average spectral difference of all image overlap areas in the initial stitched image, where the spectral difference is the difference between the spectra of each pixel point of the two sub-detection areas;

[0103] Obtain the brightness value of each pixel in the two sub-detection areas in the overlapping area through the optical satellite images of the sub-detection areas, and calculate the average brightness difference of all overlapping areas in the initial stitched image. The brightness difference is the difference in brightness of each pixel in the two sub-detection areas.

[0104] Key points are extracted from the overlapping areas of each sub-detection area using the Harris corner detection method. The key points of the two overlapping sub-detection areas are matched using the FLANN method. The Euclidean distance of the matching points is calculated using the Euclidean distance. The average Euclidean distance of all matching points in the initial stitched image is calculated to obtain the average geometric offset error.

[0105] Harris corner detection is a classic feature point detection method used to find corners in an image. This method calculates the local structure matrix of the image based on the gradient and determines whether a point is a corner by analyzing the gradient changes of pixels within a window.

[0106] FLANN is an efficient feature matching method specifically designed for nearest neighbor search in large-scale data. In remote sensing image stitching, FLANN uses an efficient search structure to match extracted key points and find corresponding points in adjacent images.

[0107] The average spectral difference, average brightness difference, and average geometric offset error are normalized, and the stitching error index is calculated based on the normalized average spectral difference, average brightness difference, and average geometric offset error. The specific acquisition steps are as follows:

[0108] SEI=ln(1+SD)+BD+e GO ;

[0109] Where SEI is the stitching error index, SD is the average spectral difference, BD is the average brightness difference, and GO is the average geometric offset error.

[0110] In this embodiment, it should be specifically explained that the steps of determining the stitching consistency of the initial stitched images according to the stitching error index are as follows:

[0111] The stitching error index is compared with the error threshold. If the stitching error index is less than the error threshold, the initial stitching image is judged to be stitched consistently; if the stitching error index is greater than or equal to the error threshold, the initial stitching image is judged to be stitched inconsistently. The error threshold is obtained by the adaptive threshold method.

[0112] Step 6: If it is determined that the initial stitched images are inconsistent, the initial stitched images are stitched and optimized to obtain the actual stitched image at the acquisition time point, thereby ensuring the overall consistency of the images and improving the accuracy of the land use classification. It should be noted that stitching optimization of the initial stitched images is a prior art, and this embodiment does not describe the specific steps in detail.

[0113] Stitching optimization is a technique used to improve the consistency of remote sensing image stitching. It aims to reduce spectral differences in overlapping areas, eliminate sudden brightness changes at the stitching boundaries, and correct geometric misalignments to improve the overall continuity and visual quality of the imagery. Stitching optimization typically includes processing methods such as color balancing, radiometric normalization, image fusion, and geometric correction.

[0114] Step 7: Store the actual stitched images at each acquisition time point into the database and perform land use classification using deep learning methods.

[0115] The deep learning method is an intelligent algorithm based on neural networks that can automatically learn features in images for high-precision land use classification. Specifically, a common method is the convolutional neural network, which can extract land feature information such as cultivated land, forests, water bodies, buildings, etc. from remote sensing images. First, the actual stitched images stored in the database are preprocessed, including radiation correction, denoising, enhancement and other operations, and then the convolutional neural network is used to extract features and divide the images into different categories. During the training phase, the model is supervised by using labeled remote sensing image data so that it can automatically identify land types. Finally, the trained deep learning model is used to perform pixel-level classification on the new image to generate a land use classification map, thereby achieving efficient and automated land use monitoring.

[0116] In this embodiment, it should be specifically explained that, Figure 2 As shown, a land use classification system based on remote sensing technology includes:

[0117] The benchmark data acquisition module divides the target area into multiple sub-detection areas, selects a benchmark time point, and collects multiple remote sensing data for each sub-detection area at the benchmark time point. The collected multiple remote sensing data are used as the benchmark data of the sub-detection area. The remote sensing data includes optical satellite images, radar images, and elevation images. The benchmark data includes benchmark optical satellite images, benchmark radar images, and benchmark elevation images. The benchmark data is then transmitted to the illumination angle change judgment module.

[0118] A multiple remote sensing data acquisition module is used to set a fixed acquisition time period. For each sub-detection area, it acquires multiple remote sensing data in real time at the current acquisition time point. The real-time multiple remote sensing data includes real-time optical satellite images and real-time radar images, and transmits the real-time multiple remote sensing data to the illumination angle change judgment module.

[0119] The illumination angle change judgment module is used to evaluate the comprehensive illumination change index of each sub-detection area based on real-time multiple remote sensing data and benchmark data, and to judge the illumination angle change based on the comprehensive illumination change index;

[0120] The illumination correction module performs illumination normalization and shadow removal on the real-time optical satellite image if it determines that the illumination angle has changed, obtaining a corrected optical satellite image. The corrected optical satellite image is then combined with the radar image and the elevation image to obtain an image of the sub-detection area, which is then transmitted to the stitching consistency judgment module.

[0121] A stitching consistency judgment module is used to stitch the images of each sub-detection area based on the same acquisition time point to obtain an initial stitching image, obtain a stitching error index based on the initial stitching image, and perform stitching consistency judgment on the initial stitching image based on the stitching error index;

[0122] The stitching optimization module performs stitching optimization on the initial stitching image to obtain the actual stitching image if it determines that the stitching of the initial stitching image is inconsistent, and transmits the actual stitching image to the land use classification module;

[0123] The land use classification module is used to store the actual mosaic images of each acquisition time point into the database and perform land use classification in combination with deep learning methods.

[0124] Finally: The above description is only a preferred embodiment of the present invention and is not intended to limit the present invention. Any modifications, equivalent substitutions, improvements, etc. made within the spirit and principles of the present invention should be included in the scope of protection of the present invention.

[0125] The above description is merely a specific embodiment of the present application, but the scope of protection of the present application is not limited thereto. Any changes or substitutions that can be easily conceived by a person skilled in the art within the technical scope disclosed in this application should be included in the scope of protection of this application. Therefore, the scope of protection of this application should be based on the scope of protection of the claims.

Claims

1. A land use classification method based on remote sensing technology, characterized in that: The following steps are involved: Step 1: Divide the target area into multiple sub-detection areas, select a benchmark time point, and collect multiple remote sensing data for each sub-detection area at the benchmark time point. The collected multiple remote sensing data are used as the benchmark data of the sub-detection area. The remote sensing data includes optical satellite images, radar images, and elevation images. The benchmark data includes benchmark optical satellite images, benchmark radar images, and benchmark elevation images. Step 2: Set a fixed acquisition time period. For each sub-detection area, obtain real-time multiple remote sensing data at the current acquisition time point. The real-time multiple remote sensing data includes real-time optical satellite images and real-time radar images. Step 3: Based on the real-time multiple remote sensing data and benchmark data, the comprehensive illumination change index of each sub-detection area is obtained, and the illumination angle change is judged based on the comprehensive illumination change index; Step 4: If it is determined that the illumination angle has changed, the real-time optical satellite image is subjected to image illumination normalization and shadow removal to obtain a corrected optical satellite image. The sub-detection area image is obtained based on the corrected optical satellite image, combined with the radar image and the elevation image; Step 5: Based on the same acquisition time point, stitch the images of each sub-detection area to obtain an initial stitched image. Evaluate the initial stitched image to obtain a stitching error index, and judge the stitching consistency of the initial stitched image based on the stitching error index. Step 6: If it is determined that the initial stitched images are not stitched consistently, then the initial stitched images are stitched and optimized to obtain the actual stitched images; Step 7: Store the actual stitched images at each acquisition time point into the database and perform land use classification using deep learning methods.

2. The land use classification method based on remote sensing technology according to claim 1, characterized in that: The steps for obtaining the comprehensive illumination change index are: Obtain the main band data of real-time optical satellite images and benchmark optical satellite images. The main band data includes the red, short-wave infrared, and near-infrared bands of each pixel point. The optical image variation coefficient is evaluated based on the main band data of the real-time optical satellite images and the benchmark optical satellite images. Obtain the polarization component of each pixel point of the real-time radar image and the reference radar image. The polarization component includes the vertical-vertical polarization component and the vertical-horizontal polarization component. The radar variation coefficient is calculated based on the polarization component of each pixel point. Obtain benchmark elevation images and evaluate the elevation influence coefficient based on the benchmark elevation images; The optical image variation coefficient, radar variation coefficient, and elevation influence coefficient are normalized, and the comprehensive illumination variation index is obtained based on the normalized optical image variation coefficient, radar variation coefficient, and elevation influence coefficient. The specific acquisition steps are as follows: Where LC is the comprehensive illumination change index, IC is the optical image change coefficient, RCC is the radar change coefficient, EIC is the elevation influence coefficient, and a1, a2, and a3 are the weight coefficients of the optical image change coefficient, the radar change coefficient, and the elevation influence coefficient.

3. The land use classification method based on remote sensing technology according to claim 2, characterized in that: The steps for obtaining the optical image variation coefficient are as follows: For each pixel, calculate the spectral change between the current time point and the reference time point. The spectral change includes the red light difference, short-wave infrared difference, and near-infrared difference. The spectral change of the pixel point is used as the feature vector, the spectral change of each pixel point is used as the data set, and the spectral change of each pixel point in the data set is used as the data point; The silhouette coefficient method is used to obtain the optimal number of clusters, and the K-means clustering method is used to cluster the data set to obtain the final clusters; The spectral variation mean of each cluster is calculated based on the final cluster, and the optical image variation coefficient is calculated based on the spectral variation mean of each cluster. The specific acquisition steps are as follows: Where IC is the optical image variation coefficient, K is the number of clusters, and C k It is expressed as the mean spectral change of the kth cluster.

4. The land use classification method based on remote sensing technology according to claim 3, characterized in that: The steps of clustering the data set using the K-means clustering method to obtain the final cluster are as follows: Step 3.1: Randomly select K data points in the data set as initial cluster centers. For each data point, calculate its Euclidean distance to each initial cluster center. For each data point, traverse the K initial cluster centers and assign it to the cluster corresponding to the initial cluster center closest to it. Step 3.2: After traversing all data points, we get the initial clusters. For each initial cluster, we calculate the mean of the data points in it and get the new cluster center. Step 3.3: Repeat steps 3.1 and 3.2 until the cluster center no longer changes, and the final cluster is obtained.

5. The land use classification method based on remote sensing technology according to claim 2, characterized in that: The radar variation coefficient acquisition step is: Obtain the vertical-vertical polarization component and vertical-horizontal polarization component of each pixel of the real-time radar image and the reference radar image; For each pixel, the logarithmic ratio change rate is calculated based on the vertical-vertical polarization component and the vertical-horizontal polarization component; The radar variation coefficient is calculated based on the logarithmic ratio change rate of the vertical-vertical polarization component and the vertical-horizontal polarization component of each pixel. The specific acquisition steps are as follows: Where RCC is the radar variation coefficient, N is the number of pixels, Δσ 0,i Expressed as the logarithmic ratio change rate of the vertical-vertical polarization component of the i-th pixel, Δσ 1,i It is expressed as the logarithmic ratio change rate of the vertical to horizontal polarization components of the i-th pixel.

6. The land use classification method based on remote sensing technology according to claim 2, characterized in that: The steps for obtaining the elevation influence coefficient are as follows: For each pixel point in the benchmark elevation image, obtain the elevation value, and calculate the slope and slope direction of each pixel point based on the elevation value of each pixel point; Obtain a reference solar altitude angle and a reference solar azimuth angle based on a reference optical satellite image, and obtain a real-time solar altitude angle and a real-time solar azimuth angle based on a real-time optical satellite image; The solar incident angle at the reference time point is calculated based on the reference solar altitude angle and the reference solar azimuth angle, and the real-time solar incident angle is calculated based on the real-time solar altitude angle and the real-time solar azimuth angle; The elevation influence coefficient is calculated based on the solar incidence angle. The specific steps are as follows: Where EIC is the elevation influence coefficient, N is the number of pixels, and θ i (T ref ) is the solar incident angle of the i-th pixel at the reference time point, θ i (T) is the real-time solar incident angle of the i-th pixel.

7. The land use classification method based on remote sensing technology according to claim 1, characterized in that: The steps of determining the change of illumination angle according to the comprehensive illumination change index are as follows: The comprehensive illumination change index is compared with the illumination change threshold. If the comprehensive illumination change index is less than the illumination change threshold, it is determined that the illumination angle has not changed; if the comprehensive illumination change index is greater than or equal to the illumination change threshold, it is determined that the illumination angle has changed.

8. The land use classification method based on remote sensing technology according to claim 1, characterized in that: The steps for obtaining the splicing error index are: The images of each sub-detection area are stitched together to obtain an initial stitched image. The image overlapping area in the initial stitched image is obtained. The spectral value of each pixel point of the two sub-detection areas in the image overlapping area is obtained through the optical satellite images of the sub-detection areas. The average spectral difference of all the image overlapping areas in the initial stitched image is calculated; Obtain the brightness value of each pixel in the two sub-detection areas in the overlapping area through the optical satellite images of the sub-detection areas, and calculate the average brightness difference of all overlapping areas in the initial stitched image; Key points are extracted from the overlapping areas of each sub-detection area using the Harris corner detection method. The key points of the two overlapping sub-detection areas are matched using the FLANN method. The Euclidean distance of the matching points is calculated using the Euclidean distance. The average Euclidean distance of all matching points in the initial stitched image is calculated to obtain the average geometric offset error. The average spectral difference, average brightness difference, and average geometric offset error are normalized, and the stitching error index is calculated based on the normalized average spectral difference, average brightness difference, and average geometric offset error. The specific acquisition steps are as follows: SIX=ln(1+SD)+BD+e GO ; Where SEI is the stitching error index, SD is the average spectral difference, BD is the average brightness difference, and GO is the average geometric offset error.

9. The land use classification method based on remote sensing technology according to claim 1, characterized in that: The step of judging the consistency of the initial stitched images according to the stitching error index is as follows: The stitching error index is compared with the error threshold. If the stitching error index is less than the error threshold, the initial stitching image is judged to be stitched consistently; if the stitching error index is greater than or equal to the error threshold, the initial stitching image is judged to be stitched inconsistently.

10. A land use classification system based on remote sensing technology, used to implement the land use classification method based on remote sensing technology according to any one of claims 1 to 9, characterized in that: The system comprises: The benchmark data acquisition module divides the target area into multiple sub-detection areas, selects a benchmark time point, and collects multiple remote sensing data for each sub-detection area at the benchmark time point. The collected multiple remote sensing data are used as the benchmark data of the sub-detection area. The remote sensing data includes optical satellite images, radar images, and elevation images. The benchmark data is then transmitted to the illumination angle change judgment module. A multiple remote sensing data acquisition module is used to set a fixed acquisition time period. For each sub-detection area, it acquires multiple remote sensing data in real time at the current acquisition time point. The real-time multiple remote sensing data includes real-time optical satellite images and real-time radar images, and transmits the real-time multiple remote sensing data to the illumination angle change judgment module. The illumination angle change judgment module is used to evaluate the comprehensive illumination change index of each sub-detection area based on real-time multiple remote sensing data and benchmark data, and to judge the illumination angle change based on the comprehensive illumination change index; The illumination correction module performs illumination normalization and shadow removal on the real-time optical satellite image if it determines that the illumination angle has changed, obtaining a corrected optical satellite image. The corrected optical satellite image is then combined with the radar image and the elevation image to obtain an image of the sub-detection area, which is then transmitted to the stitching consistency judgment module. A stitching consistency judgment module is used to stitch the images of each sub-detection area based on the same acquisition time point to obtain an initial stitching image, obtain a stitching error index based on the initial stitching image, and perform stitching consistency judgment on the initial stitching image based on the stitching error index; The stitching optimization module performs stitching optimization on the initial stitching image to obtain the actual stitching image if it determines that the stitching of the initial stitching image is inconsistent, and transmits the actual stitching image to the land use classification module; The land use classification module is used to store the actual mosaic images of each acquisition time point into the database and perform land use classification in combination with deep learning methods.

Images