A karst rocky desertification remote sensing interpretation method based on lithological background and land use preprocessing

Abstract

The present application relates to remote sensing image processing and geographic information system technical field, disclose a kind of based on lithology background and land use preprocessing Karst rocky desertification remote sensing interpretation method, comprising: S1: based on the logical constraint preprocessing of lithology background and land use, generate three types of output area, including: effective interpretation area, exclusion area, and non-karst area;S2: adaptive spatial preprocessing is carried out to effective interpretation area, using multi-scale fusion filtering and through the closed-loop adjustment driven by image quality index, output final preprocessing image;S3: in the effective interpretation area of final preprocessing image, using improved CART decision tree is classified;S4: the classified result is post-processed, and then the classified result after post-processing is spliced with the exclusion area and the non-karst area that have been assigned, to generate the final rocky desertification spatial distribution map.The present application significantly improves the accuracy and efficiency of Karst rocky desertification interpretation.

Description

Technical Field

[0001] This invention relates to the field of remote sensing image processing and geographic information system technology, specifically to a remote sensing interpretation method for karst rocky desertification based on lithological background and land use preprocessing. Background Technology

[0002] Karst desertification is a major ecological problem hindering the sustainable development of southwestern my country. Utilizing remote sensing technology for large-scale, long-term monitoring and assessment of desertification is currently the mainstream approach. Among existing technologies, the CART decision tree-based desertification interpretation method is relatively common. This method classifies the spectral and textural features of remote sensing images by constructing a decision tree model to extract different levels of desertification information.

[0003] However, existing methods have the following core flaws:

[0004] 1) Lack of preprocessing logic leads to physical paradoxes in interpretation results: Traditional methods usually classify directly on global images, failing to effectively remove "non-karst areas" and "land types that cannot be desertified in karst areas" during the preprocessing stage. This results in a large workload for interpretation and often leads to basic errors such as incorrectly labeling desertified patches on non-carbonate rock strata or identifying urban construction land as desertified areas.

[0005] 2) Sensitive to image noise and ground feature boundaries, resulting in an imbalance between interpretation accuracy and efficiency: Raw remote sensing images contain sensor noise, and the boundaries of rocky desertification patches are fragmented with blurred transition zones. Traditional CART models directly classify raw pixels, which easily misclassifies noise as information, and at the same time, it is difficult to accurately capture complex boundaries, resulting in a large tree structure and low interpretation efficiency.

[0006] 3) The splitting criteria are too simplistic and the ability to identify extreme rocky desertification areas is insufficient: Conventional CART uses overall purity (such as mean square error) as the optimization target, which is not sensitive to extremely severe rocky desertification areas with small sample sizes but serious hazards.

[0007] 4) Axial parallel splitting leads to bloated tree structure: Traditional CART uses axial parallel splitting, which requires generating extremely deep tree structures to achieve satisfactory accuracy when features are correlated or decision boundaries are oblique.

[0008] To address the aforementioned shortcomings, this invention not only introduces lithology-land use logic constraints in the preprocessing stage, but also makes systematic improvements to the core aspects of decision tree classification (splitting criteria, splitting methods, and sample weights). Both parts have equally important innovative value and together constitute the complete technical solution of this invention. Summary of the Invention

[0009] This invention aims to overcome the shortcomings of existing technologies and provide a high-precision, high-efficiency, and logically consistent remote sensing interpretation method for karst rocky desertification. This invention employs an integrated process of "adaptive spatial preprocessing with lithology-land use constraints + improved CART decision tree classification" to exclude physically impossible rocky desertification features at the source, perform targeted noise reduction and boundary enhancement on the preserved areas, and then combine this with decision tree classification that considers class imbalance and spatial consistency to achieve rapid and high-precision interpretation of rocky desertification information.

[0010] To achieve the above objectives, the following technical solution is adopted:

[0011] This invention provides a remote sensing interpretation method for karst desertification based on lithological background and land use preprocessing, comprising the following steps:

[0012] S1: Based on the logical constraints of lithological background and land use, preprocessing generates three types of output areas, including: effective interpretation areas, exclusion areas directly assigned as non-rocky desertification areas, and non-karst areas assigned as non-applicable areas.

[0013] S2: Adaptive spatial preprocessing is performed on the effective interpretation region, multi-scale fusion filtering is adopted and closed-loop adjustment driven by image quality index is used to output the final preprocessed image;

[0014] S3: Within the effective interpretation area of ​​the final preprocessed image, an improved CART decision tree is used for classification to obtain the classification result; the improved CART decision tree includes: feature space construction with physical constraints, node selection criteria guided by priority identification of high-level rocky desertification, and sample weighting mechanism;

[0015] S4: Post-process the classification results, and then stitch the post-processed classification results with the excluded areas and non-karst areas that have been assigned values ​​in step S1 to generate the final spatial distribution map of rocky desertification.

[0016] Furthermore, step S1 specifically includes:

[0017] Obtain digital geological maps of the study area to generate karst potential area masks;

[0018] Obtain land use data for the study area to extract land types that are unlikely to be affected by rocky desertification and generate exclusion masks;

[0019] The karst potential area mask and the exclusion mask are superimposed to define the three types of output regions: the effective interpretation region is the region within the karst potential area mask that does not belong to the exclusion mask; the exclusion region directly assigned as non-rocky desertification is the region corresponding to the exclusion mask; and the non-karst region assigned as non-applicable region is the region outside the karst potential area mask.

[0020] Furthermore, in step S2, the adaptive spatial preprocessing includes:

[0021] The local mean and local variance are calculated by sliding window, and the smoothing weight is determined based on the normalized local variance. The smoothing result of window mean and the enhancement result of anti-sharpening mask are weighted and combined to obtain the single-scale filter output.

[0022] The single-scale filtered output is calculated using sliding windows of different sizes, and then weighted and fused to obtain a multi-scale preprocessed image.

[0023] Calculate the average gray-level gradient of the images before and after preprocessing. When the average gray-level gradient of the multi-scale preprocessed image does not fall into the target range set based on the average gray-level gradient before preprocessing, adjust the sensitivity coefficient in the smoothing weight. After multiple consecutive adjustments of the sensitivity coefficient are ineffective, adjust the multi-scale fusion weight and iteratively generate a new multi-scale preprocessed image until the condition is met, and output the final preprocessed image.

[0024] Further, in step S2, the calculation of local mean and local variance through a sliding window, the determination of smoothing weights based on the normalized local variance, and the weighted combination of the window mean smoothing result and the anti-sharpening mask enhancement result to obtain a single-scale filter output specifically include:

[0025] Calculate the local mean and local variance of pixels within a sliding window, and use the ratio of local variance to global variance as the normalized local variance.

[0026] The smoothing weights are calculated based on the normalized local variance using a weighting function, which is: ;in, Sensitivity coefficient To normalize the local variance, This is the weighting function, used to calculate the weighting coefficients between the smoothing and sharpening components;

[0027] The mean value of the pixels within the window is calculated as the window mean smoothing result, and the difference between the original pixel value and its window mean smoothing result is calculated, multiplied by the sharpening intensity parameter, and then added to the original pixel value to obtain the anti-sharpening mask enhancement result.

[0028] The window mean smoothing result and the anti-sharpening mask enhancement result are weighted and combined according to the smoothing weight to obtain a single-scale filter output, wherein the larger the smoothing weight, the higher the proportion of the smoothing result.

[0029] Furthermore, in step S3, the feature space that integrates physical constraints includes: multi-band reflectance, normalized vegetation index, texture features, carbonate rock subclass coding, slope, and land use subclass coding that may be subject to rocky desertification.

[0030] The texture features are calculated as follows: the gray-level co-occurrence matrix is ​​calculated based on the near-infrared band, and a 7×7 sliding window, four-directional averaging and 32-level gray-level quantization are used to extract the mean and contrast as two-dimensional texture features.

[0031] Furthermore, in the node selection criterion guided by the priority identification of high-grade rocky desertification, the comprehensive scoring function of candidate splits is defined as the weighted sum of the decrease in weighted mean square error and the increase in the concentration of high-grade rocky desertification samples.

[0032] The decrease in the weighted mean square error is equal to the weighted mean square error of the parent node minus the weighted sum of the weighted mean square errors of the left and right child nodes, where the weighting coefficients of the left and right child nodes are the proportions of the sample size of the child node to the sample size of the parent node, respectively; wherein, the weighted mean square error is calculated based on the total weight of each sample within the node and its rocky desertification level, and the higher the total weight of the sample, the greater its contribution to the mean square error.

[0033] The increase in concentration of high-grade rocky desertification samples is the sum of the increase in the proportion of high-grade samples in the left child node relative to the parent node and the increase in the proportion of high-grade samples in the right child node relative to the parent node, with zero when the increment is negative.

[0034] When a candidate split results in any child node being empty, that candidate split is not included in the scoring.

[0035] Furthermore, the comprehensive scoring function for the candidate split is calculated according to the following formula:

[0036]

[0037] in, For candidate splits Overall score The decrease in the weighted mean square error. The increase in concentration of high-grade rocky desertification samples with a rocky desertification level ≥4; This is the preference coefficient, used to balance the weights of the decrease in weighted mean square error and the increase in the concentration of high-grade rocky desertification samples in the comprehensive score. This is the index of the current tree node; The current tree node The sample set; The set of left child node samples formed after the candidate split; The set of right child node samples formed after the candidate split; The number of samples in the sample set of the left child node; The number of samples in the sample set of the right child node; : The number of samples in the current node's sample set; : The set of high-level samples with a rocky desertification level ≥4 in the left child node; : The set of high-level samples with a rocky desertification level ≥4 in the right child node; : The set of high-level samples with a rocky desertification level ≥4 in the current node; The number of high-level samples in the left child node; The number of high-level samples in the right child node; The number of high-level samples in the current node.

[0038] Furthermore, in the sample weighting mechanism, the total weight of each sample is equal to the product of the base weight and the confidence weight; the base weight is taken as the desertification level of the sample, and the higher the desertification level, the greater the base weight; the confidence weight is determined according to the following objective judgment criteria: for verification points that have been measured in the field and confirmed by location, a confidence weight of 1.5 is assigned; for basic training samples, pre-classification is first performed using an unweighted traditional CART decision tree to obtain the initial classification of each sample, and then different confidence weights are assigned according to the different distance intervals that the spatial distance from the sample point to the boundary of the initial classification falls into; wherein, a value of 0.5 is assigned when the distance is less than 1 pixel, a value of 0.8 is assigned when the distance is between 1 and 3 pixels, and a value of 1.0 is assigned when the distance is greater than 3 pixels.

[0039] Furthermore, the improved CART decision tree also includes a finite form of oblique splitting mechanism, specifically including: at the current node, selecting multiple continuous features with the highest correlation to the rocky desertification level, solving the linear combination coefficients through linear discriminant analysis, and constructing synthetic features for node splitting; setting a fallback mechanism: when the sample size of the current node is less than or equal to the number of selected features plus two, or when either the high rocky desertification level group or the low rocky desertification level group is empty, it directly degenerates into axis-parallel splitting; when the condition number of the covariance matrix exceeds a preset threshold, ridge regularization is applied first, and if the condition number still exceeds the preset threshold after regularization, it degenerates into axis-parallel splitting.

[0040] Furthermore, in step S4, the post-processing includes: applying a preset minimum patch area threshold and spatial neighborhood consistency rule to the classification result to eliminate salt-and-pepper noise, and then stitching the noise-eliminated classification result with the excluded area and the non-karst area to generate the final spatial distribution map of rocky desertification.

[0041] Compared with the prior art, the present invention has the following beneficial effects:

[0042] (1) Improvement of preprocessing logic: The lithological background and land use rules are embedded as pre-processing logical constraints, which significantly reduces the misclassification phenomenon that contradicts the geological background and basically eliminates the misclassification of rocky desertification outside the carbonate rock distribution area within the study area. At the same time, a dual structure of "mask hard removal + feature soft constraint" is formed in the decision tree through lith and lucc features.

[0043] (2) Adaptation of image enhancement technology: In view of the technical characteristics of fragmented and blurred boundaries of karst rocky desertification patches, this invention introduces adaptive spatial filtering and multi-scale edge enhancement mechanism into the preprocessing stage. Through closed-loop adjustment driven by GMG relative index, the preprocessing process can be quantitatively evaluated and automatically adjusted.

[0044] (3) Targeted optimization of splitting criteria: A node scoring function Score(s) is proposed to prioritize the identification of extreme rocky desertification, so that the tree growth process changes from minimizing the overall error to optimizing the recall rate of severely rocky desertification areas. Through a sample weighting mechanism with objective criteria (pre-classification + spatial distance), the risk of misclassification or missed identification of severely rocky desertification patches is effectively reduced.

[0045] (4) Extension of splitting method: Controlled oblique splitting is introduced, and a small sample and empty group fallback mechanism is set to avoid the problem of excessive tree structure caused by traditional axial parallel splitting, while ensuring the robustness of the algorithm under extreme node conditions.

[0046] It should be understood that the description in the Summary of the Invention is not intended to limit the key or essential features of the embodiments of the present invention, nor is it intended to restrict the scope of the invention. Other features of the invention will become readily apparent from the following description. Attached Figure Description

[0047] The above and other features, advantages, and aspects of the various embodiments of the present invention will become more apparent from the accompanying drawings and the following detailed description. The drawings are provided for a better understanding of the invention and are not intended to limit the scope of the invention. In the drawings, the same or similar reference numerals denote the same or similar elements, wherein:

[0048] Figure 1This is a flowchart illustrating the steps of a remote sensing interpretation method for karst desertification based on lithological background and land use preprocessing, according to an embodiment of the present invention.

[0049] Figure 2 This is an overall technical flowchart of a remote sensing interpretation method for karst rocky desertification based on lithological background and land use preprocessing, according to an embodiment of the present invention.

[0050] Figure 3 This is a diagram comparing the effects of single-scale and multi-scale filtering;

[0051] Figure 4 This is a schematic diagram of the internal structure of the splitting node of the improved CART decision tree according to an embodiment of the present invention. Detailed Implementation

[0052] To make the objectives, technical solutions, and advantages of the embodiments of the present invention clearer, the technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0053] Furthermore, the term "and / or" in this article is merely a description of the relationship between related objects, indicating that three relationships can exist. For example, A and / or B can represent: A existing alone, A and B existing simultaneously, or B existing alone. Additionally, the character " / " in this article generally indicates that the preceding and following related objects have an "or" relationship.

[0054] Example 1:

[0055] Figure 1 This is a flowchart illustrating the steps of a remote sensing interpretation method for karst desertification based on lithological background and land use preprocessing, according to an embodiment of the present invention. Figure 2 This is a complete technical flowchart, fully illustrating the main steps and decision branches from input to output. The method includes the following steps:

[0056] S1: Based on the logical constraints of lithological background and land use, preprocessing generates three types of output areas, including: effective interpretation areas, exclusion areas directly assigned as non-rocky desertification areas, and non-karst areas assigned as non-applicable areas.

[0057] Step S1 is used to implement logical constraint preprocessing based on lithological background and land use. Specifically, it includes: acquiring a digital geological map of the study area to generate a karst potential area mask; acquiring land use data of the study area to extract land cover layers where rocky desertification is impossible, generating an exclusion mask; and overlaying the karst potential area mask with the exclusion mask to define the three types of output areas: the effective interpretation area is the area within the karst potential area mask but not belonging to the exclusion mask; the exclusion area directly assigned as non-rocky desertification is the area corresponding to the exclusion mask; and the non-karst area assigned as non-applicable area is the area outside the karst potential area mask. Specifically, it includes the following steps:

[0058] S1.1: Obtain digital geological maps of the study area, extract the distribution areas of carbonate rock strata, and generate karst potential area masks. The original remote sensing images were cut to remove non-carbonate rock areas.

[0059] S1.2: Obtain land use data for the study area, extract the "land types where rocky desertification is impossible" layer, which includes at least: large water bodies, urban construction land, road land, and perennial paddy fields, and generate an exclusion mask. .

[0060] S1.3: Superimpose the results of S1.1 and S1.2, and define three types of output regions:

[0061] ;

[0062] in: : Karst potential zone mask (carbonate rock strata distribution area); Rocky desertification cannot be covered by land types (large bodies of water, urban construction land, road land, perennial paddy fields, etc.). Effective interpretation of region masks; : The complement (i.e., the area excluding land categories); : The complement (i.e., non-carbonate rock areas); subsequent classification steps only occur in the supplementary set (i.e., non-carbonate rock areas); It will be carried out within the region.

[0063] S2: Adaptive spatial preprocessing is performed on the effective interpretation region, multi-scale fusion filtering is adopted and closed-loop adjustment driven by image quality index is used to output the final preprocessed image;

[0064] Step S2 is used to implement adaptive spatial preprocessing based on lithological constraints. The adaptive spatial preprocessing includes: calculating local mean and local variance using a sliding window, determining smoothing weights based on normalized local variance, and weighting the window mean smoothing result and the anti-sharpening mask enhancement result to obtain a single-scale filtered output; calculating the single-scale filtered output using sliding windows of different sizes and performing weighted fusion to obtain a multi-scale preprocessed image; calculating the average gray-level gradient of the images before and after preprocessing; when the average gray-level gradient of the multi-scale preprocessed image does not fall into the target range set based on the average gray-level gradient before preprocessing, adjusting the sensitivity coefficient in the smoothing weights; and after multiple ineffective adjustments to the sensitivity coefficient, adjusting the multi-scale fusion weights and iteratively generating a new multi-scale preprocessed image until the conditions are met, outputting the final preprocessed image. Specifically, it includes the following steps:

[0065] S2.1: Potential rocky desertification areas after treatment in step S1 The adaptive spatial filter is initialized using default parameters.

[0066] The process involves calculating local mean and local variance using a sliding window, determining smoothing weights based on the normalized local variance, and then weighting and combining the window mean smoothing result and the anti-sharpening mask enhancement result to obtain a single-scale filtered output. Specifically, this includes: calculating the local mean and local variance of pixels within the sliding window, and calculating the ratio of the local variance to the overall image variance as the normalized local variance; calculating the smoothing weights based on the normalized local variance using a weighting function; calculating the mean of pixels within the window as the window mean smoothing result, and multiplying the difference between the original pixel value and its window mean smoothing result by a sharpening intensity parameter, then adding this difference to the original pixel value to obtain the anti-sharpening mask enhancement result; and weighting and combining the window mean smoothing result and the anti-sharpening mask enhancement result according to the smoothing weights to obtain a single-scale filtered output. The larger the smoothing weight, the higher the proportion of the smoothed result. The detailed process is as follows:

[0067] With sliding window Calculate the local mean of pixels within the sliding window, using units of 1. and local variance . :by A sliding window centered on a pixel, with a window size of [value missing]. ,in, The side length of the sliding window (default value 7).

[0068] Filtered output Defined as a weighted combination of the smoothing and sharpening components:

[0069]

[0070] in: : Filtered output, i.e., the center pixel The grayscale value after smoothing and sharpening weighted combination;

[0071] Smoothing results for window mean: ;

[0072] Enhancement results for the anti-sharpening mask: ; The sharpening intensity parameter is set to a fixed value of 1.0.

[0073] The weighting function is used to calculate the weighting coefficients for the smoothing and sharpening components.

[0074]

[0075] in, This is the sensitivity coefficient, ranging from 0.1 to 2.0, with an initial default value of [value missing]. . This is to normalize the local variance (divided by the overall variance). : Global variance (used to normalize local variance).

[0076] For Landsat surface reflectance data (normalized to [0,1]), the default parameters are valid; for other dimensional data... It should be adjusted proportionally to the dynamic range of the image, so that The median falls within the range of 0.5 to 2. The sensitivity coefficient β is used to control the sensitivity of the filter to local variance. The larger β is, the faster the smoothing weight decays as the local variance increases, and the stronger the edge preservation ability; conversely, the smaller β is, the more the filter tends to be globally smooth.

[0077] S2.2: Introduce a multi-scale response fusion mechanism.

[0078] Use small windows respectively (Default 3x3) and large window (Default 11×11) Calculate the filtering result and The current preprocessed image is obtained through weighted fusion:

[0079]

[0080] in, Small window, default size 3×3; Large window, default size 11×11; Use a small window The calculated filtering results; Use large window The calculated filtering results; The current preprocessed image obtained after multi-scale weighted fusion; Multi-scale fusion weights (initial default value 0.5).

[0081] S2.3: Introduce an automatic evaluation and closed-loop adjustment mechanism for preprocessing quality.

[0082] Define the average gray-scale gradient (GMG) as a quality metric:

[0083]

[0084] in, Average gray-level gradient, used as an image quality evaluation indicator; The number of rows in the image; The number of columns in the image.

[0085] The absolute value depends on the image dimensions and spatial resolution. This invention uses a relative index: let the image before preprocessing be... for After pretreatment Set the target interval as The upper limit of 1.5 is based on experimental verification that when... When the increase exceeds 50%, the misclassification rate caused by false edges in subsequent classifications rises significantly; therefore, the upper limit is set to [missing value]. .

[0086] Comparative experiments conducted by the inventors involved setting filter parameters for GMG amplification of 30%, 50%, 70%, and 100% for images of the study area, and evaluating classification accuracy using the same classifier. The results showed that when the GMG amplification was within 50%, the overall classification accuracy remained above 90%; when the amplification exceeded 50%, the misclassification rate introduced by false edges increased by an average of approximately 5-8 percentage points; and when the amplification reached 100%, the misclassification rate increased by more than 15 percentage points. Therefore, the upper limit of the GMG target interval was set to 1.5 times the value before preprocessing.

[0087] The closed-loop adjustment process is as follows:

[0088] (1) Calculate the current of value;

[0089] (2) If Accept the present As the final preprocessed image Exit S2;

[0090] (3) If (Overly smoothed) preferentially increase (Step 0.1, maximum 2.0), maintain the currently adjusted... If the value remains unchanged, after 3 consecutive adjustments If it does not recover, it will begin to enlarge. (Step 0.05, maximum 0.7), re-execute S2.1-S2.2 to generate new... ;

[0091] (4) If (Over-sharpening), prioritize reducing (Step 0.1, lower limit 0.1), maintain the currently adjusted... If the value remains unchanged, after 3 consecutive adjustments If it still hasn't recovered, it will begin to decrease. (Step 0.05, lower limit 0.3), re-execute S2.1-S2.2 to generate new... .

[0092] (5) Iterate until If the interval requirement is met or the maximum number of iterations (default 10) is reached, and if the requirement is still not met after 10 iterations, then... Output the set of parameters that is closest to the midpoint of the interval.

[0093] in, : The average grayscale gradient value of the image before preprocessing; The final output preprocessed image (meets the requirements) (Target range requirements).

[0094] The selection of the above step size and limit value is based on the dynamic range (0-1) of the Landsat series satellite surface reflectance products and the image characteristics of typical karst landform areas, which are determined through grid search optimization. For other sensors or image types, those skilled in the art can adjust them accordingly based on the dynamic range of the image and the complexity of the ground features.

[0095] Figure 3 This is a schematic diagram comparing the effects of single-scale and multi-scale filtering. Figure 3 The comparison of surface reflectance curves along the horizontal profile (y=105) of the simulated karst landform image is shown: the original image (blue dotted line) shows obvious high-frequency noise fluctuations; the single-scale adaptive filtering (red dotted line) has a certain degree of noise suppression but tends to be overly smooth; the multi-scale fusion filtering (green solid line) effectively suppresses noise while better preserving the reflectance gradient jump at the bare rock-vegetation boundary, and the edge details are clearer.

[0096] S3: Within the effective interpretation area of ​​the final preprocessed image, an improved CART decision tree is used for classification to obtain the classification result; the improved CART decision tree includes: feature space construction with physical constraints, node selection criteria guided by priority identification of high-level rocky desertification, and sample weighting mechanism;

[0097] Step S3: Improved CART decision tree classification

[0098] S3.1: Construct a feature space that incorporates physical constraints.

[0099] The feature space incorporating physical constraints includes: multi-band reflectance, normalized vegetation index, texture features, carbonate rock subclass coding, slope, and land use subclass coding for land prone to desertification. The texture features are calculated as follows: based on the near-infrared band, a gray-level co-occurrence matrix is ​​calculated, using a 7×7 sliding window, four-directional averaging, and 32-level gray-level quantization; the mean and contrast are then extracted as two-dimensional texture features. Specifically:

[0100] for Extract the feature vector for each pixel within the region. :

[0101]

[0102] in, For reflectivity in each band, Normalized Difference Vegetation Index (NDVI), texture is the textural feature, lith is the code for carbonate rock subclass (such as limestone, dolomite, carbonate rock interbedded with clastic rocks, etc.), slope is the slope, and lucc is the code for land use subclass that may be subject to rocky desertification (such as dry land, shrubland, bare land, sparse grassland, etc.).

[0103] The specific calculation specifications for the texture feature are as follows: The gray-level co-occurrence matrix (GLCM) is calculated based on the near-infrared band (which is most sensitive to differences between rocks and vegetation), using a 7×7 sliding window and four directions (…). The average and 32-level grayscale were used to extract two statistical measures, the mean and contrast, as two-dimensional texture features (denoted as ). and To enhance the characterization of details in rocky desertification patches.

[0104] S3.2: Introduce a node selection criterion guided by the priority identification of high-level rocky desertification.

[0105] In the node selection criterion guided by the priority identification of high-grade rocky desertification, the comprehensive scoring function for candidate splits is defined as the weighted sum of the decrease in weighted mean square error and the increase in the concentration of high-grade rocky desertification samples. The decrease in weighted mean square error is equal to the weighted mean square error of the parent node minus the weighted mean square errors of the left and right child nodes, where the weighting coefficients of the left and right child nodes are the proportions of the child node's sample size to the parent node's sample size, respectively. The weighted mean square error is calculated based on the total weight of each sample within the node and its rocky desertification level; the higher the total weight of the sample, the greater its contribution to the mean square error. The increase in the concentration of high-grade rocky desertification samples is the sum of the increment of the proportion of high-grade samples in the left child node relative to the parent node and the increment of the proportion of high-grade samples in the right child node relative to the parent node, with zero when the increment is negative. When a candidate split results in any child node being empty, that candidate split does not participate in the scoring. Specifically:

[0106] Set the current node The sample set is Rocky desertification level Among them: 1 = no rocky desertification, 2 = potential rocky desertification, 3 = slight rocky desertification, 4 = moderate rocky desertification, 5 = extremely severe rocky desertification. Candidate splitting Divide the node into its left child node. and right child node When a candidate split results in any child node being empty ( or When this happens, the candidate split is skipped and does not participate. Calculation. The comprehensive scoring function is defined as follows:

[0107]

[0108] in, Candidate splitting schemes Overall score; The preference coefficient is used to balance the weights of the decrease in weighted mean square error and the increase in the concentration of high-grade rocky desertification samples in the comprehensive score. In this embodiment, it is taken as... (Parameter sensitivity analysis shows) (The effect is stable at the time).

[0109] The decrease in weighted mean square error:

[0110]

[0111] in, : The index of the current tree node; : Current tree node The sample set; : Candidate splitting schemes; The set of left child node samples formed after the candidate split; The set of right child node samples formed after the candidate split; : The number of samples in the sample set of the left child node; : The number of samples in the sample set of the right child node; : The number of samples in the current node's sample set; Candidate splitting schemes The resulting decrease in weighted mean square error; Sample set The weighted mean square error.

[0112] To increase the concentration of high-grade rocky desertification samples (grade ≥ 4):

[0113]

[0114] Candidate splitting schemes The resulting increase in the concentration of high-grade rocky desertification samples; : The set of high-level samples with a rocky desertification level ≥4 in the left child node; : The set of high-level samples with a rocky desertification level ≥4 in the right child node; : The set of high-level samples with a rocky desertification level ≥4 in the current node; The number of high-level samples in the left child node; The number of high-level samples in the right child node; The number of high-level samples in the current node.

[0115] S3.3: Introduce sample weighting and cost-sensitive mechanisms.

[0116] In the sample weighting mechanism, the total weight of each sample is equal to the product of the base weight and the confidence weight. The base weight is taken as the desertification level of the sample; the higher the desertification level, the greater the base weight. The confidence weight is determined according to the following objective judgment criteria: for verification points that have been measured in the field and confirmed by location, a confidence weight of 1.5 is assigned. For basic training samples, pre-classification is first performed using an unweighted traditional CART decision tree to obtain the initial classification of each sample. Then, based on the different distance intervals that the sample point falls into from the boundary of the initial classification, different confidence weights are assigned. The specific process is as follows:

[0117] For each sample Assign weights :

[0118]

[0119] in, : No. Feature vectors of each sample; : No. The rocky desertification level of each sample; : No. The total weight of each sample; The basic weighting function takes values ​​representing the sample's desertification level. (The higher the level, the greater the weight); : No. Field / boundary confidence weights for each sample (determined based on field validation or distance to class boundaries).

[0120] in, (The higher the level of rocky desertification, the greater the weight). The objective judgment criteria are:

[0121] Verification points measured in the field and confirmed by GPS positioning ;

[0122] For the basic training samples, firstly, an unweighted traditional CART decision tree is used to... Perform a pre-classification of the region to obtain the initial class of each sample; calculate the spatial distance (in pixels) from the sample point to the boundary of the initial class: if the distance is less than 1 pixel, assign... When the distance is between 1 and 3 pixels, assign When the distance is greater than 3 pixels, assign ;

[0123] Other general samples are taken .

[0124] The formula for calculating weighted node purity (MSE) is:

[0125]

[0126] Sample set The weighted average level of rocky desertification.

[0127] S3.4: Introduce a finite form of oblique splitting mechanism and set a small sample fallback.

[0128] Finite form of oblique splitting mechanism: At the current node, select multiple continuous features with the highest correlation to the rocky desertification level, solve the linear combination coefficient through linear discriminant analysis, and construct synthetic features for node splitting; set a fallback mechanism: when the sample size of the current node is less than or equal to the number of selected features plus two, or when either the high rocky desertification level group or the low rocky desertification level group is empty, directly degenerate to axis-parallel splitting; when the condition number of the covariance matrix exceeds a preset threshold, ridge regularization is applied first. If the condition number still exceeds the preset threshold after regularization, then degenerate to axis-parallel splitting. The specific process is as follows:

[0129] When splitting a node, first calculate the full sample size of the current node. Each feature and The correlations are analyzed, and the top correlations are selected. The continuous features (including reflectance, NDVI, texture features (including mean and contrast), and slope, excluding discrete coded features lith and lucc) are defined. (Default value is 3). Select the linear combination coefficient vector. Constructing synthetic features:

[0130]

[0131] : Linear combination coefficient vector ; : No. The linear combination coefficients of the features; : No. The first sample One eigenvalue; :Depend on A synthetic feature constructed by linear combination of features.

[0132] The method for classifying high and low rocky desertification levels is as follows: calculate the current node. The median of all rocky desertification grades in the sample; grades greater than The samples were classified as high-level rocky desertification group ( ), the level is less than or equal to The samples were classified as low-level rocky desertification group ( This transforms the multi-class classification problem into a binary linear discriminant analysis at the current node. The calculation of high and low rocky desertification levels is performed in the same... The feature mean vector in the feature space With covariance matrix .coefficient The objective function is obtained by solving the weighted linear discriminant analysis:

[0133]

[0134] : Current node The median of the rocky desertification level across all samples; High-level rocky desertification group, consisting of levels greater than […]. The composition of the sample; Low-level rocky desertification group, consisting of levels less than or equal to The composition of the sample; High rocky desertification level group in the selected The feature mean vector in each feature space; Low-level rocky desertification group in the selected The feature mean vector in each feature space; : Covariance matrix of high rocky desertification grade groups; : Covariance matrix of low-level rocky desertification groups.

[0135] Fallback mechanism:

[0136] If the current node's sample size Or high rocky desertification level group With low rocky desertification level group When any set in the set is empty, it directly degenerates into a traditional axial parallel split (K=1) without attempting regularization.

[0137] If the condition number of the covariance matrix exceeds When applying ridge regularization, priority is given to applying it. Then recalculate the condition number; if the condition number still exceeds the limit after regularization... If so, it degenerates into axial parallel splitting (K=1). Ridge regularization parameter, with values ​​ranging from 1 to 2. It is used to improve the invertibility of the covariance matrix.

[0138] The node information records the dominant characteristics and their coefficients, which can be used by technicians to trace the basis for splitting. When K=1, oblique splitting degenerates into axial parallel splitting, which is fully compatible with traditional CART.

[0139] Unlike general oblique decision trees, the oblique splitting mechanism of this invention is specially designed for the special characteristics of karst rocky desertification classification: ① oblique splitting is performed only on continuous features, while discrete features remain on parallel axes to avoid overfitting; ② the number of features participating in oblique splitting is controlled by correlation screening to reduce computational complexity; ③ a multi-layer backoff mechanism is set to ensure that the algorithm can still run stably under extremely small sample or ill-conditioned matrix conditions.

[0140] S3.5: Recursively execute S3.2 to S3.4 to construct the complete CART decision tree. A post-pruning strategy based on the weighted bias reduction rate is adopted: when... When the division stops, The default value is 0.05 (parameter sensitivity analysis shows...). The pruning effect is stable at [0.03, 0.08]. : Current node The collective term for the set of all child node samples after splitting (i.e.) and (the union of the two sets of the same set). The threshold for the reduction rate of deviation after pruning is preset to 0.05 (parameter sensitivity analysis shows that the pruning effect is stable when the value range is 0.03~0.08).

[0141] S4: Post-process the classification results, and then stitch the post-processed classification results with the excluded areas and non-karst areas that have been assigned values ​​in step S1 to generate the final spatial distribution map of rocky desertification.

[0142] Step S4 is used to perform post-processing and output results.

[0143] Combining the three output regions from step S1 and the classification results from step S3, a minimum patch area threshold is applied. (For example, 9 pixels, i.e., 3×3) and spatial neighborhood consistency rules (majority filtering) are used to eliminate salt-and-pepper noise. Then, the classification result after noise removal is compared with the exclusion region already assigned in step S1. Non-karst areas The images are stitched together to generate the final spatial distribution map of rocky desertification, without assigning duplicate values.

[0144] Example 2:

[0145] The technical solution of the present invention will be described in detail below with reference to an example from the Yinjiang River Basin (karst trough area) in Guizhou Province.

[0146] Example: Extraction of information on rocky desertification in the Yinjiang River Basin, Guizhou Province

[0147] 1. Data Preparation: Landsat remote sensing images of the study area from multiple historical periods were collected, covering both TM and OLI sensors. The surface reflectance product (SR-ARD) provided by USGS was used uniformly, and radiometric normalization was performed across different image periods. GMG calculations directly used the raw reflectance dimensions (range 0–1, dimensionless). Simultaneously, 1:50,000 digital geological maps, contemporaneous land use data, and SRTM DEM (30m resolution) were collected. Field validation data included 320 field survey points, of which 45 were severely rocky desertification sites, 68 were severely rocky desertification sites, 87 were slightly rocky desertification sites, 65 were potential rocky desertification sites, and 55 were non-rocky desertification sites. The time difference between the survey and image acquisition did not exceed 6 months.

[0148] 2. Execute step S1 (logical constraint preprocessing): Generate a karst area mask based on the geological map. Extracting water bodies, towns, and paddy fields from land use data to generate exclusion masks. .definition The area (excluding the carbonate rock area) will proceed to subsequent processing; The region is directly assigned the value of "non-rocky desertification (human-made / water body)"; the region of non-carbonate rock is assigned the value of "non-applicable area".

[0149] 3. Perform step S2 (adaptive spatial preprocessing):

[0150] Initial settings: sliding window Sharpening intensity .

[0151] Before computational preprocessing of images The reflectivity range is 0~1. Taking a certain period of Landsat 8 OLI near-infrared band as an example, 0.223. The Yinjiang River basin is a karst peak-cluster depression landform with dramatic topographic relief and a strong mixture of bare rock and vegetation, resulting in significant spectral differences between adjacent pixels. This value is within a reasonable range.

[0152] generate ,calculate Falling within the target range No parameter adjustments are needed; output is performed directly. .

[0153] If GMG does not meet the interval requirements, adjust sequentially according to the step rule in S2.3. and Iterate until the target is met.

[0154] 4. Perform step S3 (Improved CART classification):

[0155] Sample preparation: Based on visual interpretation of high-resolution images and field validation, a total of 1500 training samples were collected, including 1180 basic training samples (obtained from high-resolution image interpretation) and 320 field survey points (not overlapping with the basic samples). The field points were also used for subsequent independent validation.

[0156] Pre-classification and sample weight calculation: First, an unweighted traditional CART decision tree is used to classify the samples. A coarse classification is performed on the region to obtain the initial class of each basic training sample. The spatial distance from each basic sample point to the boundary of the initial class is calculated: among the 1180 basic samples, 23 samples have a distance <1 pixel (assigned a distance of 1 pixel to each basic sample point). ), 67 samples of pixels 1-3 (assigned) ), 1090 samples with a resolution > 3 pixels (assigned) ) Weights were directly assigned to 320 field survey sampling points. .

[0157] Feature construction: extract Landsat 1-5 band reflectance, NDVI, texture features (calculated based on near-infrared band GLCM mean and contrast, 2 dimensions in total), lithological subclass coding (limestone, dolomite, carbonate rock with clastic rocks 3 types), slope, and land use subclass coding (dryland, shrubland, bare land, sparse grassland 4 types).

[0158] Weighted training: setting sample weights Preference coefficient Deviation reduction rate threshold Oblique split K=2 (only continuous features: reflectivity, NDVI, texture, slope) and processed according to the fallback priority of step S3.4.

[0159] The final decision tree contains 47 nodes and has a depth of 9.

[0160] 5. Perform step S4 (post-processing and output): Apply the minimum patch area threshold. (3×3 pixels) and 3×3 majority filtering are used to output the final rocky desertification level distribution map.

[0161] IV. Ablation Experiment and Effect Verification

[0162] Experimental hardware environment: Intel Core i7-10750H CPU@2.60GHz, 16GB RAM; Software: R language 4.0.5 + rpart package. The processing time statistics range is from reading the preprocessed image (S1 mask generation has been completed; S1 is offline preprocessing, and the results can be reused in multiple periods) to outputting the classification result map, including feature extraction, model training, and full-map prediction, excluding the time spent on S1 GIS mask generation (this step is completed offline and does not affect the online processing time comparison). The baseline method (traditional CART) uses rpart's default parameters, only uses spectral features, and does not include the S1 mask, S2 filtering, or S3 improvement.

[0163] The comparison results of the various schemes are shown in Table 1:

[0164] Table 1

[0165] Each module contributes independently to the overall accuracy. The complete combination of the three modules achieves an overall accuracy of 91.3% and a recall rate of 86.5% for extremely high recall, both superior to any single or dual-module combination. The complete method saves only 3 minutes compared to the +S1+S2 scheme. This is mainly because the pre-classification step in the S3 improvement adds an extra full-map prediction time (approximately 2 minutes), which is partially offset by the approximately 5-minute saving brought about by the significant reduction in tree depth (from 19 layers to 9 layers), resulting in a net saving of 3 minutes. The number of misclassifications contradicting the geological background is reduced from 27 to 3, essentially eliminating misclassifications of rocky desertification outside carbonate rock distribution areas.

[0166] The modules of this invention can be selectively activated according to the actual application scenario: for images with high data quality and clear ground feature boundaries, the S1+S3 scheme can be used alone; for images with severe noise, the S2 module can be activated; for scenarios requiring high-precision identification of extremely severe rocky desertification, the complete S1+S2+S3 scheme can be used. The modular design of this invention makes it highly adaptable to different scenarios and is not a redundant design.

[0167] Parameter sensitivity description: Cross-validation testing was performed. Between 0.2 and 0.4 Between 0.03 and 0.08, Within the range of 0.3 to 0.6, the overall accuracy fluctuation does not exceed ±1.5%, indicating that the method is not sensitive to parameters. The best results are achieved when K=2, K=1 increases the tree depth by about 30%, and K=3 increases the risk of overfitting.

[0168] This invention discloses a remote sensing interpretation method for karst rocky desertification based on lithological background and land use preprocessing. The method first generates a logical constraint mask using lithological distribution and land use data, limiting the interpretation scope to the effective area within the carbonate rock region that does not exclude land types. Then, adaptive multi-scale filtering is applied to the effective area for image enhancement, and parameters are automatically optimized through closed-loop adjustment driven by gray-level average gradient. Finally, in an improved CART decision tree, multi-source physical features are integrated, a splitting criterion that considers high-level rocky desertification identification is adopted, a sample weighting mechanism is introduced, and finite-form oblique splitting is used to complete high-precision classification. This invention effectively solves the problems of physical paradoxes, noise sensitivity, insufficient identification of extreme rocky desertification, and bloated tree structures in traditional methods, significantly improving the accuracy and efficiency of rocky desertification interpretation.

[0169] The above description is only a preferred embodiment of the present invention and is used only to explain the present invention. It 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 within the protection scope of the present invention.

[0170] It should also be noted that, in the embodiments of this application, relational terms such as "first" and "second" are used only to distinguish one entity or operation from another entity or operation, and do not necessarily require or imply any such actual relationship or order between these entities or operations. Furthermore, the terms "comprising," "including," or any other variations thereof are intended to cover non-exclusive inclusion, such that a process, method, article, or apparatus that comprises a list of elements includes not only those elements but also other elements not expressly listed, or elements inherent to such a process, method, article, or apparatus. Without further limitations, an element defined by the phrase "comprising one..." does not exclude the presence of other identical elements in the process, method, article, or apparatus that includes said element.

[0171] The above description of the disclosed embodiments enables those skilled in the art to make or use this application. Various modifications to these embodiments will be readily apparent to those skilled in the art, and the general principles defined in the embodiments of this application may be implemented in other embodiments without departing from the spirit or scope of this application. Therefore, this application is not to be limited to the embodiments shown in this application, but is to be accorded the widest scope consistent with the principles and novel features disclosed in the embodiments of this application.

Claims

1. A remote sensing interpretation method for karst desertification based on lithological background and land use preprocessing, characterized in that, Includes the following steps: S1: Based on the logical constraints of lithological background and land use, preprocessing generates three types of output areas, including: effective interpretation areas, exclusion areas directly assigned as non-rocky desertification areas, and non-karst areas assigned as non-applicable areas. S2: Adaptive spatial preprocessing is performed on the effective interpretation region, multi-scale fusion filtering is adopted and closed-loop adjustment driven by image quality index is used to output the final preprocessed image; S3: Within the effective interpretation area of ​​the final preprocessed image, an improved CART decision tree is used for classification to obtain the classification result; the improved CART decision tree includes: feature space construction with physical constraints, node selection criteria guided by priority identification of high-level rocky desertification, and sample weighting mechanism; S4: Post-process the classification results, and then stitch the post-processed classification results with the excluded areas and non-karst areas that have been assigned values ​​in step S1 to generate the final spatial distribution map of rocky desertification.

2. The method according to claim 1, characterized in that, Step S1 specifically includes: Obtain digital geological maps of the study area to generate karst potential area masks; Obtain land use data for the study area to extract land types that are unlikely to be affected by rocky desertification and generate exclusion masks; The karst potential area mask and the exclusion mask are superimposed to define the three types of output regions: the effective interpretation region is the region within the karst potential area mask that does not belong to the exclusion mask; the exclusion region directly assigned as non-rocky desertification is the region corresponding to the exclusion mask; and the non-karst region assigned as non-applicable region is the region outside the karst potential area mask.

3. The method according to claim 1, characterized in that, In step S2, the adaptive spatial preprocessing includes: The local mean and local variance are calculated by sliding window, and the smoothing weight is determined based on the normalized local variance. The smoothing result of window mean and the enhancement result of anti-sharpening mask are weighted and combined to obtain the single-scale filter output. The single-scale filtered output is calculated using sliding windows of different sizes, and then weighted and fused to obtain a multi-scale preprocessed image. Calculate the average gray-level gradient of the images before and after preprocessing. When the average gray-level gradient of the multi-scale preprocessed image does not fall into the target range set based on the average gray-level gradient before preprocessing, adjust the sensitivity coefficient in the smoothing weight. After multiple consecutive adjustments of the sensitivity coefficient are ineffective, adjust the multi-scale fusion weight and iteratively generate a new multi-scale preprocessed image until the condition is met, and output the final preprocessed image.

4. The method according to claim 3, characterized in that, In step S2, the calculation of local mean and local variance using a sliding window, the determination of smoothing weights based on normalized local variance, and the weighted combination of the window mean smoothing result and the anti-sharpening mask enhancement result to obtain a single-scale filter output specifically include: Calculate the local mean and local variance of pixels within a sliding window, and use the ratio of local variance to global variance as the normalized local variance. The smoothing weights are calculated based on the normalized local variance using a weighting function, which is: ;in, Sensitivity coefficient To normalize the local variance, This is the weighting function, used to calculate the weighting coefficients between the smoothing and sharpening components; The mean value of the pixels within the window is calculated as the window mean smoothing result, and the difference between the original pixel value and its window mean smoothing result is calculated, multiplied by the sharpening intensity parameter, and then added to the original pixel value to obtain the anti-sharpening mask enhancement result. The window mean smoothing result and the anti-sharpening mask enhancement result are weighted and combined according to the smoothing weight to obtain a single-scale filter output, wherein the larger the smoothing weight, the higher the proportion of the smoothing result.

5. The method according to claim 1, characterized in that, In step S3, the feature space that integrates physical constraints includes: multi-band reflectance, normalized vegetation index, texture features, carbonate rock subclass coding, slope, and land use subclass coding that may be subject to rocky desertification. The texture features are calculated as follows: the gray-level co-occurrence matrix is ​​calculated based on the near-infrared band, and the mean and contrast are extracted as two-dimensional texture features using a 7×7 sliding window, four-directional averaging and 32-level gray-level quantization.

6. The method according to claim 1, characterized in that, In the node selection criterion guided by the priority identification of high-grade rocky desertification, the comprehensive scoring function of candidate splits is defined as the weighted sum of the decrease in weighted mean square error and the increase in the concentration of high-grade rocky desertification samples. The decrease in the weighted mean square error is equal to the weighted mean square error of the parent node minus the weighted sum of the weighted mean square errors of the left and right child nodes, where the weighting coefficients of the left and right child nodes are the proportions of the sample size of the child node to the sample size of the parent node, respectively; wherein, the weighted mean square error is calculated based on the total weight of each sample within the node and its rocky desertification level, and the higher the total weight of the sample, the greater its contribution to the mean square error. The increase in concentration of high-grade rocky desertification samples is the sum of the increase in the proportion of high-grade samples in the left child node relative to the parent node and the increase in the proportion of high-grade samples in the right child node relative to the parent node, with zero when the increment is negative. When a candidate split results in any child node being empty, that candidate split is not included in the scoring.

7. The method according to claim 6, characterized in that, The comprehensive scoring function for the candidate splits is calculated according to the following formula: ； in, For candidate splits Overall score The decrease in the weighted mean square error. The increase in concentration of high-grade rocky desertification samples with a rocky desertification level ≥4; This is the preference coefficient, used to balance the weights of the decrease in weighted mean square error and the increase in the concentration of high-grade rocky desertification samples in the comprehensive score. This is the index of the current tree node; The current tree node The sample set; The set of left child node samples formed after the candidate split; The set of right child node samples formed after the candidate split; : The number of samples in the sample set of the left child node; : The number of samples in the sample set of the right child node; : The number of samples in the current node's sample set; : The set of high-level samples with a rocky desertification level ≥4 in the left child node; : The set of high-level samples with a rocky desertification level ≥4 in the right child node; : The set of high-level samples with a rocky desertification level ≥4 in the current node; The number of high-level samples in the left child node; The number of high-level samples in the right child node; The number of high-level samples in the current node.

8. The method according to claim 1, characterized in that, In the sample weighting mechanism, the total weight of each sample is equal to the product of the base weight and the confidence weight; The basic weight is taken as the rocky desertification level of the sample, and the higher the rocky desertification level, the greater the basic weight. The confidence weights are determined based on the following objective criteria: For verification points that have been measured in the field and confirmed by location, a confidence weight of 1.5 is assigned; For basic training samples, pre-classification is first performed using an unweighted traditional CART decision tree to obtain the initial classification of each sample, and then different confidence weights are assigned according to the different distance intervals that the sample point falls into from the boundary of the initial classification; where the value is 0.5 when the distance is less than 1 pixel, 0.8 when the distance is between 1 and 3 pixels, and 1.0 when the distance is greater than 3 pixels.

9. The method according to claim 1, characterized in that, The improved CART decision tree also includes a finite number of oblique splitting mechanisms, specifically: At the current node, select several continuous features that are most correlated with the level of rocky desertification, solve the linear combination coefficients through linear discriminant analysis, and construct synthetic features for node splitting; Set a fallback mechanism: when the sample size of the current node is less than or equal to the number of selected features plus two, or when either the high rocky desertification level group or the low rocky desertification level group is empty, it will directly degenerate into axial parallel splitting; When the condition number of the covariance matrix exceeds When applying ridge regularization, the regularization parameter of the ridge regularization is given priority. If the condition number still exceeds the limit after regularization Then it degenerates into axial parallel splitting.

10. The method according to claim 1, characterized in that, In step S4, the post-processing includes: applying a preset minimum patch area threshold and spatial neighborhood consistency rule to the classification result to eliminate salt-and-pepper noise, and then stitching the noise-eliminated classification result with the excluded area and the non-karst area to generate the final spatial distribution map of rocky desertification.

Images