Land utilization classification method and system based on remote sensing technology

Abstract

The invention relates to the field of land classification, and discloses a land utilization classification method and system based on a remote sensing technology, which are used for solving the problem that the illumination angle is changed and is not matched with the boundary of a splicing region during land utilization classification, and comprises the following steps: collecting reference data of a sub-detection region; acquiring various real-time remote sensing data for each sub-detection area, evaluating to obtain a comprehensive illumination change index, performing illumination angle change judgment, if the illumination angle is judged to change, performing image illumination normalization and shadow removal to obtain sub-detection area images, and performing image splicing on each sub-detection area image to obtain a sub-detection area image; the method comprises the following steps: obtaining an initial spliced image, evaluating to obtain a splicing error index, performing splicing consistency judgment according to the splicing error index, and if the splicing of the initial spliced image is judged to be inconsistent, performing splicing optimization on the initial spliced image to obtain an actual spliced image, thereby effectively improving the accuracy of remote sensing image acquisition and improving the rationality of land utilization classification.

Description

Technical Field

[0001] The present invention relates to the field of land classification, and more specifically 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 important tool for analyzing land cover types and changes. Through remote sensing images, researchers and decision makers can efficiently perform tasks such as land classification. Land use classification based on remote sensing images usually involves the analysis of multi-temporal data in order to identify dynamic changes in land types and provide accurate classification results.

[0003] Existing remote sensing land use classification methods usually rely on the acquisition of multi-source data, including optical satellite images, radar images, and elevation images, and are 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, the time interval for acquiring remote sensing images is usually fixed. 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 solution 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: In practical applications, even if images are collected at fixed time intervals, remote sensing images at different times may still show inconsistent representations of the same area in different images due to factors such as lighting conditions and changes in observation angles. For example, due to changes in the solar altitude angle, the image brightness and shadow length of the same object at different times may change, affecting the stability of the classification algorithm. In addition, due to slight offsets in the satellite orbit or different imaging angles, the boundaries of the stitched areas may not match. This error not only affects the consistency of land use classification results, but may also lead to misjudgment of surface cover types, thus affecting subsequent analysis and decision-making. Summary of the invention

[0005] 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.

[0006] To achieve the above object, the present invention provides the following technical solutions: A land use classification method based on remote sensing technology includes the following steps: Step 1: Divide a 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 in combination with the radar image and the elevation image; Step 5: Based on the same acquisition time point, each sub-detection area image is stitched to obtain an initial stitched image, and a stitching error index is obtained based on the initial stitched image evaluation, and the initial stitched image is judged for stitching consistency based on the stitching error index; Step 6: If it is judged that the initial stitched image is inconsistent, the initial stitched image is optimized for stitching to obtain an actual stitched image; Step 7: The actual stitched images of each acquisition time point are stored in the database, and land use classification is performed in combination with the deep learning method.

[0007] Preferably, the step of acquiring the comprehensive illumination change index is: acquiring main band data of the real-time optical satellite image and the reference optical satellite image, the main band data including the red light, short-wave infrared and near-infrared bands of each pixel, evaluating the optical image change coefficient according to the main band data of the real-time optical satellite image and the reference optical satellite image; acquiring 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, calculating the radar change coefficient according to the polarization component of each pixel; acquiring the reference elevation image, evaluating the elevation influence coefficient according to the reference elevation image; normalizing the optical image change coefficient, the radar change coefficient and the elevation influence coefficient, evaluating the comprehensive illumination change index according to the normalized optical image change coefficient, the radar change coefficient and the elevation influence coefficient, the specific acquisition steps are: ; In the formula, Expressed as a comprehensive light change index, Expressed as the optical image variation coefficient, Expressed as the radar variation coefficient, Expressed as the elevation influence coefficient, It is the weight coefficient of the optical image variation coefficient, the weight coefficient of the radar variation coefficient and the weight coefficient of the elevation influence coefficient.

[0008] Preferably, the optical image variation coefficient acquisition step is: for each pixel point, calculate the spectral variation between the current time point and the reference time point, the spectral variation includes red light difference, short-wave infrared difference and near infrared difference; take the spectral variation of the pixel point as a feature vector, take the spectral variation of each pixel point as a data set, and the spectral variation of each pixel point in the data set as a data point; use the silhouette coefficient method to obtain the optimal number of clusters, use the K-means clustering method to cluster the data set to obtain the final cluster cluster; calculate the spectral variation mean of each cluster cluster according to the final cluster cluster, and calculate the optical image variation coefficient according to the spectral variation mean of each cluster cluster. The specific acquisition steps are: ; In the formula, It is expressed as the optical image variation coefficient, K is the number of clusters, It is expressed as the mean spectral change of the kth cluster.

[0009] Preferably, the steps of clustering the data set using the K-means clustering method to obtain the final cluster clusters 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 therein 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.

[0010] Preferably, the radar variation coefficient acquisition step is: acquiring the vertical-vertical polarization component and the vertical-horizontal polarization component of each pixel point of the real-time radar image and the reference radar image; For each pixel point, the logarithmic ratio change rate is calculated according to the vertical-vertical polarization component and the vertical-horizontal polarization component respectively; the radar variation coefficient is calculated according to the logarithmic ratio change rate of the vertical-vertical polarization component and the vertical-horizontal polarization component of each pixel point. The specific acquisition steps are as follows: ; In the formula, It is expressed as the radar variation coefficient, N is the number of pixels, It is expressed as the logarithmic ratio change rate of the vertical-vertical polarization component of the i-th pixel, It is expressed as the logarithmic ratio change rate of the vertical-horizontal polarization component of the i-th pixel.

[0011] Preferably, the step of acquiring the elevation influence coefficient is: for each pixel point in the reference elevation image, acquiring the elevation value, and calculating the slope and slope direction of each pixel point according to the elevation value of each pixel point; acquiring the reference solar altitude angle and the reference solar azimuth angle according to the reference optical satellite image, and acquiring the real-time solar altitude angle and the real-time solar azimuth angle according to the real-time optical satellite image; calculating the solar incident angle at the reference time point according to the reference solar altitude angle and the reference solar azimuth angle, and calculating the real-time solar incident angle according to the real-time solar altitude angle and the real-time solar azimuth angle; calculating the elevation influence coefficient according to the solar incident angle, and the specific acquisition steps are: ; In the formula, It is expressed as the elevation influence coefficient, N is the number of pixels, is the solar incident angle of the i-th pixel at the reference time point, is the real-time solar incident angle of the i-th pixel.

[0012] Preferably, the step of judging the change of illumination angle according to the comprehensive illumination change index is as follows: comparing the comprehensive illumination change index with the illumination change threshold; if the comprehensive illumination change index is less than the illumination change threshold, then judging that the illumination angle has not changed; if the comprehensive illumination change index is greater than or equal to the illumination change threshold, then judging that the illumination angle has changed.

[0013] Preferably, the stitching error index acquisition step is: stitching the images of each sub-detection area to obtain an initial stitching image, obtaining the image overlap area in the initial stitching image, respectively 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 stitching image; respectively 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 stitching image; extracting key points in the image overlap area of ​​each sub-detection area through the Harris corner point detection method, matching the key points of the two overlapping sub-detection areas through the FLANN method, using the Euclidean distance to calculate the Euclidean distance of the matching points, and averaging the Euclidean distances of all matching points in the initial stitching image to obtain the average geometric offset error; normalizing the average spectral difference, the average brightness difference and the average geometric offset error, and calculating the stitching error index according to the normalized average spectral difference, the average brightness difference and the average geometric offset error. The specific acquisition steps are: ; In the formula, Expressed as the splicing error index, Expressed as the average spectral difference, Expressed as the average brightness difference.

[0014] 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, judging that the initial stitched image is stitched consistently; if the stitching error index is greater than or equal to the error threshold, judging that the initial stitched image is stitched inconsistently.

[0015] Preferably, a land use classification system based on remote sensing technology comprises: 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, 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, which is used to set a fixed collection time period, and for each sub-detection area, obtains 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, and transmits the real-time multiple remote sensing data to the illumination angle change judgment module; an illumination angle change judgment module, which 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 performs 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, combined 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 the stitching error index is obtained according to the initial stitched image evaluation, and the initial stitched image is judged for stitching consistency according to the stitching error index; the stitching optimization module, if it is judged that the initial stitched image is inconsistent, the initial stitched image is optimized for stitching 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.

[0016] Technical effects and advantages of the present invention: Collect benchmark data of sub-detection areas, obtain real-time multiple remote sensing data for each sub-detection area, evaluate and obtain the 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 image of the sub-detection area, stitch each sub-detection area image to obtain the initial stitched image, evaluate and obtain the stitching error index, and judge the stitching consistency according to 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, which effectively improves the accuracy of remote sensing image acquisition and improves the rationality of land use classification. BRIEF DESCRIPTION OF THE DRAWINGS

[0017] Figure 1 A flow chart of a land use classification method based on remote sensing technology provided in the embodiment of the present application 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

[0018] The technical solution of the present invention will be clearly and completely described below in conjunction with the accompanying drawings in the present invention. In addition, the forms of the various structures recorded in the following embodiments are only 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 recorded in the following embodiments. All other implementations obtained by ordinary technicians in this field without making creative work belong to the scope of protection of the present invention.

[0019] 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: 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 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; When the target area is divided into a plurality of sub-detection areas, a division method may be selected according to actual conditions, and the division methods include division based on fixed grids, division based on administrative areas, and division based on natural geographical units; Fixed grid division is divided according to fixed longitude and latitude grids. Common grid sizes are , as well as ; Administrative divisions are based on administrative units (such as provinces, cities, counties or townships); natural geographical units are based on natural geographical features such as watershed distribution, mountains, plains and wetlands.

[0020] In remote sensing image analysis, the selection of 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 the time point with the best image quality, stable lighting conditions, and minimal cloud interference.

[0021] In this embodiment, it should be specifically explained that the steps of selecting the reference time point are: Determine the key factors that affect image quality, including lighting conditions (moderate solar altitude angle), meteorological conditions (no cloud or low cloud cover), ground object stability (such as vegetation growth, urban expansion), and the availability of remote sensing data (optical satellite images, radar images, and elevation images are available at the same time). Select the time point that meets all the key factors as the initial reference time point; Obtain multi-temporal remote sensing data of the target area in the past 2 to 3 years. Multi-temporal remote sensing data includes optical satellite images, radar images, elevation images, and meteorological data, which are used to analyze changes in light, cloud cover, and surface features at different time points; Clouds can seriously affect the availability of remote sensing images, so it is necessary to first screen images with less cloud cover. The scene classification map of optical satellite images can be used to automatically detect cloud cover and remove images with cloud cover exceeding 10%. At the same time, candidate images are manually checked to ensure that the images are not affected by thin clouds, haze or atmospheric scattering to ensure image clarity and availability; The solar altitude directly affects the brightness and shadow range of the image, so it is necessary to select the time point with a moderate solar altitude. The solar angle information can be extracted through the remote sensing image metadata, and images with a solar altitude between 40°-70° can be selected to reduce the impact of long shadows and light overexposure. At the same time, a curve of the solar altitude change throughout the year is drawn to ensure that the lighting conditions at the selected time point are stable; The stability of surface features is crucial for illumination normalization and subsequent classification, so key indicators such as vegetation index, building index, and water index need to be calculated to screen the time point with the smallest surface change. For example, the stability of vegetation growth is calculated through the vegetation index, and the time point with the smallest vegetation index change is selected to reduce the spectral change error caused by factors such as crop growth and forest seasonal changes; Check whether there is a corresponding radar image at the initial reference time point for shadow detection and illumination 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; Taking all the screening factors into consideration, a time point that meets the following conditions is selected from all the initial benchmark time points: no cloud or low cloud cover, moderate solar altitude angle, stable surface reflectivity, and optical images, radar images, and elevation images are all available. Finally, the image quality is manually reviewed to ensure that the benchmark image has no human interference (such as road construction, pollution, etc.). This time point will be used as the benchmark time point for the sub-detection area and will be used for subsequent image stitching, illumination normalization, and land use classification analysis.

[0022] Step 2: Set a fixed collection time period, and for each sub-detection area, obtain real-time multiple remote sensing data at the current collection time point, and the real-time multiple remote sensing data include real-time optical satellite images and real-time radar images; Step 3: Obtain the comprehensive illumination change index of each sub-detection area based on the real-time multiple remote sensing data and benchmark data, and determine the illumination angle change based on the comprehensive illumination change index; In this embodiment, it should be specifically explained 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 include red light, short-wave infrared and near-infrared bands of each pixel point. The optical image variation coefficient is obtained based on the main band data of real-time optical satellite images and benchmark optical satellite images. The polarization components of each pixel of the real-time radar image and the reference radar image are obtained. The polarization components include vertical-vertical polarization components and vertical-horizontal polarization components. The radar change coefficient is calculated based on the polarization components 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 assist in determining which changes are caused by illumination and which are caused by changes in ground objects. Obtain a benchmark elevation image, and evaluate the elevation influence coefficient based on the benchmark elevation image; 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 impact 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 become darker or form shadows. By calculating the slope and aspect using 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.

[0023] The optical image variation coefficient, radar variation coefficient and elevation influence coefficient are normalized, and the comprehensive illumination variation index is obtained by evaluating the normalized optical image variation coefficient, radar variation coefficient and elevation influence coefficient. The specific acquisition steps are as follows: ; In the formula, Expressed as a comprehensive light change index, It is expressed as the optical image variation coefficient. When the optical image variation coefficient is large, it means that the image's brightness, color, and shadow characteristics have changed significantly, which may be caused by changes in the sun's angle, weather conditions, or changes in surface reflectivity. This change will directly affect the comprehensive illumination variation index, causing significant differences in the illumination conditions of the image at different time points, thereby affecting the accuracy of illumination normalization and land use classification of remote sensing images. It is expressed as the radar variation coefficient. Since radar images are not affected by light, their changes mainly come from changes in surface humidity, roughness or structure of objects. If the radar variation coefficient is high, it means that the change of objects is the main factor, not the change of light. In this case, the comprehensive light variation index may be low. If the radar variation coefficient is low, it means that the objects are relatively stable and the change of light becomes the main influencing factor. In this case, the comprehensive light variation index may be high. It is expressed as the elevation influence coefficient. Since the solar altitude angle and solar azimuth angle will change at different time points, the area with large terrain undulation will cause significant changes in the light distribution, thus affecting the image brightness, shadow position and other features. When the elevation influence coefficient is high, it means that the change of slope and slope direction has a greater impact on the solar incident angle, causing the lighting conditions to change significantly, and the comprehensive lighting change index also increases accordingly; on the contrary, when the terrain is relatively flat, the lighting conditions are relatively stable, and the degree of lighting change in the image is small. is the weight coefficient of the optical image variation coefficient, the weight coefficient of the radar variation coefficient, and the weight coefficient of the elevation influence coefficient, and , Obtained through AHP, e.g. It can be 0.4, 0.3, 0.3.

[0024] The analytic hierarchy process is a decision analysis method that decomposes complex problems into multiple factors by establishing a hierarchical structure, and calculates weights by comparing them pairwise 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, and experts or data analysis methods can be used to assign values; then, the relative weights of each factor are calculated through consistency checks to ensure that the weights are reasonable and reliable.

[0025] The polarization component is the polarization characteristic of the electromagnetic wave when the radar wave returns to the radar receiver after interacting with the surface, mainly including vertical-vertical polarization and vertical-horizontal polarization. Vertical-vertical polarization means that the electromagnetic waves emitted and received by the radar are both vertically polarized, and is often used to detect smooth surfaces, such as water bodies and urban buildings; vertical-horizontal polarization means that the electromagnetic waves emitted by the radar are vertically polarized, but are horizontally polarized when received, and are suitable for identifying complex surface features such as vegetation and forests. Analyzing the changes in these polarization components can help distinguish changes in illumination, changes in surface humidity, or changes in the objects themselves, and improve the classification accuracy of remote sensing images.

[0026] In this embodiment, it should be specifically explained that the steps of obtaining the optical image variation coefficient are: For each pixel, calculate the spectral change between the current time point and the reference time point, the spectral change includes red light difference, short-wave infrared difference and near infrared difference; The spectral change of the pixel point is taken as the feature vector, the spectral change of each pixel point is taken as the data set, and the spectral change of each pixel point in the data set is taken 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 according to the final cluster, and the optical image variation coefficient is calculated according to the spectral variation mean of each cluster. The specific acquisition steps are as follows: ; In the formula, It is expressed as the optical image variation coefficient, K is the number of clusters, It is expressed as the mean spectral change of the kth cluster.

[0027] The silhouette coefficient method is a method used to evaluate clustering effects and determine the optimal number of clusters. It determines the rationality of clustering by measuring the compactness of data points within their clusters and the degree of separation from the nearest cluster. This method can effectively avoid 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.

[0028] 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: Step 3.1: Randomly select K data points in the data set as the initial cluster centers. For each data point, calculate its Euclidean distance to each initial cluster center. The specific method is: ,in It is expressed as the Euclidean distance from the data point to the cluster center, where Represented as data points, Represented as 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; Step 3.2: After traversing all data points, the initial clusters are obtained. For each initial cluster, the mean of the data points in it is calculated to obtain a 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.

[0029] In this embodiment, it should be specifically explained that the steps of obtaining the radar variation coefficient are: Obtain the vertical-vertical polarization component and the 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 specific acquisition steps are: ; ; In the formula, Expressed as the rate of change of the logarithmic ratio of the vertical-vertical polarization component, Expressed as the rate of change of the logarithmic ratio of the vertical-horizontal polarization component, Expressed as the vertical-vertical polarization component of the real-time radar image, Expressed as the vertical-vertical polarization component of the reference radar image, Expressed as the vertical-horizontal polarization component of the real-time radar image, Expressed as the vertical-horizontal polarization component of the reference radar image; 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 point. The specific acquisition steps are as follows: ; In the formula, Expressed as radar variation coefficient, if Smaller, e.g. , indicating that most of the changes come from illumination rather than ground changes. If it is larger, it means that the ground object itself has changed and the light has little effect. N is the number of pixels. It is expressed as the logarithmic ratio change rate of the vertical-vertical polarization component of the i-th pixel, It is expressed as the logarithmic ratio change rate of the vertical-horizontal polarization component of the i-th pixel.

[0030] In this embodiment, it should be specifically explained that the steps for obtaining the elevation influence coefficient are: For each pixel point in the benchmark elevation image, the elevation value is obtained, and the slope and slope direction of each pixel point are calculated according to the elevation value of each pixel point; Obtain the benchmark solar altitude angle and the benchmark solar azimuth angle based on the benchmark optical satellite image, and obtain the real-time solar altitude angle and the real-time solar azimuth angle based on the real-time optical satellite image. The solar altitude angle 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 angle is the angle of the sun relative to the true north direction, which determines the direction of illumination and the direction of shadow projection. 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 specific acquisition steps are: ; ; In the formula, is the solar incident angle of the i-th pixel at the reference time point, is the real-time solar incident angle of the i-th pixel, and are the reference solar altitude angle and the reference solar azimuth angle, respectively. and 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; The elevation influence coefficient is calculated according to the solar incidence angle. The specific acquisition steps are as follows: ; In the formula, It is expressed as the elevation influence coefficient, where N is the number of pixels. The difference in solar incidence angles for all pixels is calculated and the average is taken to measure the impact of terrain on illumination changes.

[0031] In this embodiment, it should be specifically explained that the steps of calculating the slope and slope direction of each pixel point according to the elevation value of each pixel point are: Get the elevation values ​​of the east-west adjacent pixels of the pixel point, and calculate the east-west slope of the pixel point. The specific steps are as follows: ; In the formula, It is expressed as the east-west slope, d is the pixel resolution, and They are respectively represented as the elevation values ​​of the east and west adjacent pixels of the pixel; Get the elevation values ​​of the north and south adjacent pixels of the pixel point, and calculate the north-south slope of the pixel point. The specific steps are as follows: ; In the formula, It is expressed as the north-south slope, d is the pixel resolution, and They are respectively expressed as the elevation values ​​of the north and south adjacent pixels of the pixel; The slope is calculated based on the east-west slope and the east-west slope. The specific steps are as follows: ; Where S is the slope, It is expressed as the east-west slope, It is expressed as north-south slope; 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: ; Where A is the slope aspect, which indicates the direction of the slope relative to the north direction, ranging from 0° to 360°, and is used to determine the impact of the sunlight incident angle.

[0032] In this embodiment, it should be specifically explained 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 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. The illumination change threshold is obtained by the adaptive threshold method. The adaptive threshold method is an algorithm that dynamically determines the threshold based on data distribution. It can automatically adapt to changes in different scenes and make the threshold more accurate and stable. When calculating the illumination change threshold, the adaptive threshold method automatically selects the optimal threshold based on the data distribution of the comprehensive illumination change index.

[0033] 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. It should be specifically noted that image illumination normalization and shadow removal are prior arts, and the specific steps are not described in detail in this embodiment. 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 deviation caused by illumination changes and terrain shadows, so that images have consistent lighting conditions at different times and in different areas, and improve the accuracy of image analysis and classification.

[0034] Light 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 lighting 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 error caused by uneven lighting.

[0035] In this embodiment, it should be specifically explained that the steps of obtaining the sub-detection area image according to the corrected optical satellite image, combined with the radar image and the elevation image are as follows: The corrected optical satellite images, radar images and elevation images are aligned through geometric correction; the radar images are combined to enhance the ground object information, and the elevation images are combined to extract the terrain information. The relationship between the optical satellite images, radar images and elevation images is learned through neural networks to generate multi-source fusion images.

[0036] 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; In this embodiment, it should be specifically explained that the steps for obtaining the splicing error index are: Perform image stitching on the images of each sub-detection area to obtain an initial stitched image, obtain an image overlap area in the initial stitched image, the image overlap area is the overlapped part 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, the spectral difference is the difference of the spectrum of each pixel point of the two sub-detection areas; The brightness value of each pixel in the two sub-detection areas in the image overlap area is obtained through the optical satellite images of the sub-detection areas, and the average brightness difference of all image overlap areas in the initial stitched image is calculated. The brightness difference is the difference in brightness of each pixel in the two sub-detection areas; The Harris corner detection method is used to extract key points in the overlapping area of ​​each sub-detection area. 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 Euclidean distances of all matching points in the initial stitched image are averaged to obtain the average geometric offset error. 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 the point is a corner by analyzing the gradient changes of the pixels in the window.

[0037] FLANN is an efficient feature matching method specially designed for nearest neighbor search in large-scale data. In remote sensing image stitching, FLANN matches the extracted key points through an efficient search structure to find the corresponding points of adjacent images.

[0038] 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: ; In the formula, Expressed as the splicing error index, Expressed as the average spectral difference, Expressed as the average brightness difference.

[0039] In this embodiment, it should be specifically explained that the steps of judging the stitching consistency of the initial stitched images according to the stitching error index are 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. The error threshold is obtained by the adaptive threshold method.

[0040] Step 6: If it is determined that the initial stitched images are not stitched consistently, the initial stitched images are stitched consistently to obtain actual stitched images 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 specifically noted that stitching consistently the initial stitched images is a prior art, and this embodiment does not describe the specific steps in detail. Stitching optimization is a technique used to improve the consistency of remote sensing image stitching, which aims to reduce the spectral differences in overlapping areas, eliminate brightness mutations at the stitching boundaries, and correct geometric misalignment to improve the overall continuity and visual quality of the image. Stitching optimization usually includes color balancing, radiometric normalization, image fusion, geometric correction and other processing methods.

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

[0042] 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, the common method is the convolutional neural network, which can extract ground feature information from remote sensing images, such as cultivated land, forests, water bodies, buildings, etc. First, the actual spliced ​​images stored in the database are preprocessed, including radiation correction, denoising, enhancement and other operations, and then the convolutional neural network is used for feature extraction to divide the images into different categories. In 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 classify the new images at the pixel level to generate a land use classification map, realizing efficient and automated land use monitoring.

[0043] In this embodiment, it should be specifically explained that Figure 2 As shown, a land use classification system based on remote sensing technology includes: The benchmark data acquisition module 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. The benchmark data includes benchmark optical satellite images, benchmark radar images, and benchmark elevation images. The benchmark data is 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, multiple real-time remote sensing data is acquired 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 the real-time multiple remote sensing data is transmitted 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 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 is determined that the illumination angle has changed, and obtains 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, and the sub-detection area image is transmitted to the splicing 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 evaluation, 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 is determined that the initial stitching image is not stitched consistently, and transmits the actual stitching image to the land use classification module; The land use classification module is used to store the actual stitched images of each acquisition time point into the database and perform land use classification in combination with deep learning methods.

[0044] 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 protection scope of the present invention.

[0045] The above is only a specific implementation of the present application, but the protection scope of the present application is not limited thereto. Any person skilled in the art who is familiar with the present technical field can easily think of changes or substitutions within the technical scope disclosed in the present application, which should be included in the protection scope of the present application. Therefore, the protection scope of the present application should be based on the protection scope 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 reference time point, collect multiple remote sensing data for each sub-detection area at the reference time point, and use the collected multiple remote sensing data as the reference data of the sub-detection area. The remote sensing data includes optical satellite images, radar images, and elevation images, and the reference data includes reference optical satellite images, reference radar images, and reference elevation images. Step 2: Set a fixed collection time period, and for each sub-detection area, obtain real-time multiple remote sensing data at the current collection time point, and the real-time multiple remote sensing data include real-time optical satellite images and real-time radar images; Step 3: Obtain the comprehensive illumination change index of each sub-detection area based on the real-time multiple remote sensing data and benchmark data, and determine the illumination angle change 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 stitching images are not consistent, then the initial stitching images are stitched and optimized to obtain the actual stitching images; Step 7: Store the actual stitched images at each acquisition time point into the database and perform land use classification in combination with 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 include red light, short-wave infrared and near-infrared bands of each pixel point. The optical image variation coefficient is obtained based on the main band data of real-time optical satellite images and benchmark optical satellite images. The polarization components of each pixel point of the real-time radar image and the reference radar image are obtained. The polarization components include vertical-vertical polarization components and vertical-horizontal polarization components. The radar variation coefficient is calculated according to the polarization components of each pixel point. Obtain a benchmark elevation image, and evaluate the elevation influence coefficient based on the benchmark elevation image; The optical image variation coefficient, radar variation coefficient and elevation influence coefficient are normalized, and the comprehensive illumination variation index is obtained by evaluating the normalized optical image variation coefficient, radar variation coefficient and elevation influence coefficient. The specific acquisition steps are as follows: ； In the formula, Expressed as a comprehensive light change index, Expressed as the optical image variation coefficient, Expressed as the radar variation coefficient, Expressed as the elevation influence coefficient, It is the weight coefficient of the optical image variation coefficient, the weight coefficient of the radar variation coefficient and the weight coefficient of the elevation influence coefficient.

3. A land use classification method based on remote sensing technology according to claim 2, characterized in that: The steps of 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 red light difference, short-wave infrared difference and near infrared difference; The spectral change of the pixel point is taken as the feature vector, the spectral change of each pixel point is taken as the data set, and the spectral change of each pixel point in the data set is taken 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 according to the final cluster, and the optical image variation coefficient is calculated according to the spectral variation mean of each cluster. The specific acquisition steps are as follows: ； In the formula, It is expressed as the optical image variation coefficient, K is the number of clusters, 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 is 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 the initial cluster centers. For each data point, calculate the Euclidean distance from 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: After traversing all data points, the initial clusters are obtained. For each initial cluster, the data points in it are calculated to 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 the 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 point. The specific acquisition steps are as follows: ； In the formula, It is expressed as the radar variation coefficient, N is the number of pixels, It is expressed as the logarithmic ratio change rate of the vertical-vertical polarization component of the i-th pixel, It is expressed as the logarithmic ratio change rate of the vertical-horizontal polarization component 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, the elevation value is obtained, and the slope and slope direction of each pixel point are calculated according to the elevation value of each pixel point; Obtain a reference solar altitude angle and a reference solar azimuth angle based on reference optical satellite images, and obtain a real-time solar altitude angle and a real-time solar azimuth angle based on real-time optical satellite images; 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 according to the solar incidence angle. The specific acquisition steps are as follows: ； In the formula, It is expressed as the elevation influence coefficient, N is the number of pixels, is the solar incident angle of the i-th pixel at the reference time point, 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 step of determining the change of illumination angle according to the comprehensive illumination change index is 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 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.

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, and 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, and the average spectral difference of all image overlapping areas in the initial stitched image is calculated; The brightness value of each pixel in the two sub-detection areas in the image overlap area is obtained through the optical satellite images of the sub-detection areas, and the average brightness difference of all the image overlap areas in the initial stitched image is calculated; The Harris corner detection method is used to extract key points in the overlapping area of ​​each sub-detection area. 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 Euclidean distances of all matching points in the initial stitched image are averaged 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: ； In the formula, Expressed as the splicing error index, Expressed as the average spectral difference, Expressed as the average brightness difference.

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 a land use classification method based on remote sensing technology as claimed in 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, collects multiple remote sensing data for each sub-detection area at the benchmark time point, uses the collected multiple remote sensing data as the benchmark data of the sub-detection area, and 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. For each sub-detection area, multiple real-time remote sensing data is acquired 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 the real-time multiple remote sensing data is transmitted 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 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 is determined that the illumination angle has changed, and obtains 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, and the sub-detection area image is transmitted to the splicing 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 evaluation, 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 is determined that the initial stitching image is not stitched consistently, and transmits the actual stitching image to the land use classification module; The land use classification module is used to store the actual stitched images of each acquisition time point into the database and perform land use classification in combination with deep learning methods.

Images