A method for optimizing the classification results of a rubber plantation based on satellite embedded (SE) features

A rubber plantation classification method was constructed by using satellite embedded features. The classification results were optimized by using dynamic thresholds and neighborhood constraints, which solved the problems of misjudgment and omission by traditional remote sensing features and achieved higher accuracy in rubber plantation classification.

CN122135225APending Publication Date: 2026-06-02RUBBER RES INST CHINESE ACADEMY OF TROPICAL AGRI SCI

Patent Information

Authority / Receiving Office
CN · China
Patent Type
Applications(China)
Current Assignee / Owner
RUBBER RES INST CHINESE ACADEMY OF TROPICAL AGRI SCI
Filing Date
2026-02-09
Publication Date
2026-06-02

AI Technical Summary

Technical Problem

In existing technologies, traditional remote sensing features have low dimensionality, which makes rubber plantations easily confused with the spectral features of other vegetation, resulting in many misclassifications, poor adaptability, difficulty in balancing the elimination of misclassified pixels and the retention of true pixels, and the initial classification results have omissions and errors.

Method used

By employing satellite embedded (SE) features, a rubber plantation feature template is constructed, pixel similarity is calculated, low-confidence pixels are removed by combining dynamic thresholding, and gap filling is performed based on neighborhood density and spectral similarity to optimize the classification results.

Benefits of technology

It effectively eliminates classification noise, recovers incorrectly removed pixels, improves the spatial continuity and spectral consistency of rubber plantation classification, and enhances classification accuracy and adaptability.

✦ Generated by Eureka AI based on patent content.

Smart Images

  • Figure CN122135225A_ABST
    Figure CN122135225A_ABST
Patent Text Reader

Abstract

This invention relates to the technical field of remote sensing image classification and vegetation monitoring, and particularly to an optimization method for rubber plantation classification results based on satellite embedded (SE) features. The method includes: acquiring satellite embedded images containing N-dimensional embedding vectors; constructing a normalized feature template vector based on rubber plantation reference samples, and calculating the cosine similarity between all pixels in the region and the template to generate a similarity layer; identifying and removing low-confidence pixels in the initial classification results using a strategy based on neighborhood density and dynamic percentile thresholds; and then, based on dual constraints of neighborhood coverage and average spectral similarity, directionally filling gaps in the refined results to generate the final optimized classification result. By utilizing high-dimensional embedded features and through data-driven dynamic thresholds and spatial-spectral joint constraints, this method effectively solves the problems of misclassification, missed classification, and poor adaptability in traditional methods, significantly improving the accuracy and spatial continuity of rubber plantation classification.
Need to check novelty before this filing date? Find Prior Art

Description

Technical Field

[0001] This invention relates to the technical field of remote sensing image classification and vegetation monitoring, and in particular to an optimization method for rubber plantation classification results based on satellite embedded (SE) features. Background Technology

[0002] As important economic forest land, the accurate extraction of the spatial distribution and area of ​​rubber plantations is of great significance for agricultural resource management, ecological environment assessment, and industrial planning. Currently, rubber plantation classification methods based on remote sensing imagery are widely used, mainly relying on traditional remote sensing features such as spectral and textural features to achieve pixel-level classification.

[0003] However, existing technologies have the following shortcomings: First, traditional remote sensing features have low dimensionality, and the spectral features of rubber plantations are easily confused with those of other vegetation (such as woodlands and orchards), resulting in many misclassified pixels in the classification results; Second, post-classification processing often uses fixed thresholds or simple neighborhood analysis, which cannot adapt to the differences in rubber plantation features in different regions and years, and it is difficult to balance "removing misclassified pixels" and "preserving true pixels"; Third, there are often missed pixels in the initial classification results, and existing filling methods lack feature constraints, which easily lead to incorrect filling and damage the accuracy and spatial continuity of the classification results.

[0004] Satellite embedded features (SE features), as a type of high-dimensional remote sensing feature, can capture deep spatial and spectral correlation information of pixels, providing new technical support for the accurate differentiation of vegetation types. Therefore, there is an urgent need for a method to optimize rubber plantation classification results by fusing satellite embedded features. Summary of the Invention

[0005] The purpose of this invention is to provide an optimization method for rubber plantation classification results based on satellite embedded (SE) features, in order to solve the problems of misjudgment, omission, and poor adaptability of traditional remote sensing features in the prior art mentioned in the background.

[0006] To achieve the above objectives, the present invention provides the following technical solution: an optimization method for rubber plantation classification results based on satellite embedded (SE) features, comprising the following steps:

[0007] S1. Satellite Embedded Image Data Acquisition and Processing: Acquire satellite embedded images covering the target area, wherein each pixel of the satellite embedded image contains an N-dimensional embedding vector;

[0008] S2. Construction of embedded feature template for rubber plantation: Based on the known rubber plantation reference sample area, the embedding vector of all pixels in the sample area is extracted, the average value is calculated and normalized to obtain the feature template vector of the rubber plantation.

[0009] S3. Calculate the similarity of pixels in the entire region: Traverse each pixel in the target region, calculate the cosine similarity between its embedding vector and the feature template vector, and generate a similarity layer for the entire region.

[0010] S4. Low-confidence pixel removal based on dynamic threshold: Based on the initial rubber plantation classification results, combined with the spatial neighborhood density of the pixels and their similarity values ​​in the similarity layer, a dynamically calculated percentile threshold is used to identify and remove low-confidence pixels in the classification results to obtain a refined classification image.

[0011] S5. Gap filling based on neighborhood density and spectral similarity: Identify gap pixels in the refined classification image, fill gap pixels that simultaneously meet the conditions of the number of classified pixels in the neighborhood, neighborhood coverage, and neighborhood average spectral similarity threshold, and generate the final optimized rubber plantation classification result.

[0012] Preferably, the specific process of constructing the rubber plantation embedded feature template in step S2 includes:

[0013] S21. Generate a binary sample mask image based on the rubber plantation reference sample area;

[0014] S22. Overlay the satellite embedded image with the sample mask image and extract the embedding vector of all pixels in the sample area;

[0015] S23. Calculate the average value of all pixel embedding vectors in the sample region to obtain the initial feature template vector;

[0016] S24. Perform L2 normalization on the initial feature template vector to obtain the feature template vector in unit vector form.

[0017] Preferably, in step S3, the cosine similarity algorithm is used to calculate the similarity value between the embedding vector of each pixel and the feature template vector. The calculation formula is as follows:

[0018]

[0019] in, For pixels The similarity value ranges from [0,1]. A higher value indicates that the pixel is more similar to the rubber plantation sample in terms of features. is the feature dimension of the satellite embedded image (SE), with a maximum value of 64; j is the dimension index; For the first The first pixel Dimensional embedding feature values.

[0020] Preferably, the specific process of low-confidence pixel removal based on dynamic threshold in step S4 includes:

[0021] S41. Convert the initial rubber plantation classification results into binary images;

[0022] S42. Perform neighborhood analysis on the binary image and count the number of pixels that have been classified as rubber plantations within the predetermined neighborhood of each pixel;

[0023] S43. Construct a set of effective pixels, wherein the effective pixels are pixels that have been classified as rubber plantations and whose number of neighbors is greater than or equal to the minimum number of neighbors (removed);

[0024] S44. Calculate the p-th percentile of the similarity value corresponding to the effective pixel set in the similarity layer, and use it as the dynamic removal threshold θ;

[0025] S45. For each pixel classified as a rubber plantation, if its number of neighbors is less than the minimum number of neighbors (removed), it is retained; if its number of neighbors is greater than or equal to the minimum number of neighbors (removed) and its similarity value is greater than or equal to the dynamic removal threshold θ, it is retained; otherwise, it is removed from the classification results.

[0026] Preferably, the specific process of filling the gaps in step S5 includes:

[0027] S51. Compare the refined classification image with the initial classification result and identify pixels that were not classified as rubber plantations as candidate gap-filling pixels.

[0028] S52. For each of the candidate filling gap pixels, the following judgment is made:

[0029] Calculate the number of pixels within its predetermined filling neighborhood that have been classified as rubber plantations. ;

[0030] Calculate the neighborhood coverage rate, i.e. The ratio of the number of possible neighbors to the number of neighbors within the filled neighborhood;

[0031] Calculate the average embedding similarity value between the gap cell and all its neighboring cells that have been classified as rubber plantations. ;

[0032] S53. If the candidate filling gap pixel simultaneously meets the following conditions, it shall be classified as a rubber plantation;

[0033] (a) ≥ Minimum neighborhood number (filled);

[0034] (b) The neighborhood coverage rate is ≥ the coverage threshold;

[0035] (c) ≥ Similarity threshold.

[0036] Preferably, the percentile p in step S44 ranges from 1 to 5.

[0037] Preferably, the range of the neighborhood to be filled in step S52 is greater than the range of the neighborhood analysis in the removal phase in step S42.

[0038] Preferably, after step S53, the method further includes the step of: limiting the filling result using an additional geographic constraint mask, wherein the geographic constraint mask includes at least one of a slope mask, an elevation mask, or a water body mask.

[0039] Compared with the prior art, the beneficial effects of the present invention are:

[0040] 1. A two-stage correction strategy of "removal first, filling later" is adopted. The first stage generates a similarity layer based on sample point features, and uses dynamic percentile threshold to identify and remove low-confidence pixels, effectively eliminating classification noise. The second stage uses dual constraints of neighborhood density and spectral similarity to perform targeted filling of the missed areas and restore the effective pixels that were wrongly removed. This design avoids the problem of poor adaptability of the fixed threshold method to different regions, and prevents overexpansion during the filling process through dual constraints.

[0041] 2. Fully utilize the high-dimensional feature representation capability of pre-trained embedding vectors, eliminating the need for manual design of classification features. Employ data-driven dynamic thresholds to improve the generalization ability of the method. The spatial-spectral joint constraint mechanism ensures the spatial continuity and spectral consistency of the correction results. Attached Figure Description

[0042] Figure 1 This is a schematic diagram of the method for correcting rubber plantation classification results using embedded satellite data according to the present invention. Detailed Implementation

[0043] The technical solutions in the embodiments of the present invention will be clearly and completely described below. 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.

[0044] This invention provides a technical solution: an optimization method for rubber plantation classification results based on satellite embedded (SE) features, comprising the following steps:

[0045] Step S1: Acquisition of satellite embedded image data

[0046] 1. Based on the geographical extent of the target area and the target year, filter the image data covering the area from the satellite embedded image dataset. The satellite embedded image (GOOGLE / SATELLITE_EMBEDDING / V1 / ANNUAL) is a multi-dimensional feature image preprocessed by a deep learning model. Each pixel contains an N-dimensional embedding vector (typical value N=64), which encodes the pixel's spectral features, texture features, and temporal variation features.

[0047] 2. Spatial mosaicking is performed on the selected images to generate annual embedded images covering the entire target area. A pixel-level synthesis strategy is used during the mosaicking process to ensure that each location obtains a valid embedding vector value.

[0048] 3. Obtain reference sample areas for the target features (such as rubber plantations). The reference samples are training sample point or area vector data of known categories, which are used for the construction of subsequent template vectors.

[0049] Step S2: Construction of embedded feature template for rubber plantation

[0050] 1. Sample Mask Generation: Based on the rubber plantation reference sample area (training sample point or area vector data of known categories), a binary mask image is generated, in which pixels within the sample area are assigned a value of 1, and pixels outside the area are assigned a value of 0.

[0051] 2. Sample feature extraction: The embedded image is superimposed on the sample mask to extract the embedding vector of all pixels in the sample area.

[0052] 3. Mean Template Calculation: The mean is used to perform regional statistics on the masked feature layer, calculating the mean of each feature band within the reference sample region to obtain a 64-dimensional average embedding vector, which serves as the initial feature template vector for the rubber plantation. The calculation formula is:

[0053]

[0054] in, This represents the total number of pixels within the sample area. For the first A 64-dimensional embedding vector for each pixel.

[0055] 4. Template Normalization: The initial template vector is L2 normalized to convert it into a unit vector, eliminating the influence of feature scale differences on similarity calculation and ensuring numerical stability. The normalization formula is:

[0056]

[0057] In the formula, The L2 norm (Euclidean norm) of the template vector T is used to ultimately generate a rubber plantation feature template in unit vector form.

[0058] Step S3: Calculate the similarity of pixels across the entire region

[0059] 1. Pixel Feature Extraction: Traverse each pixel within the target region and extract its 64-dimensional embedding vector. Since the vectors in the embedded image have been pre-normalized to unit vectors, no additional processing is required.

[0060] 2. Cosine Similarity Calculation: The cosine similarity algorithm is used to calculate the similarity between each pixel's embedding vector and the rubber plantation feature template; for unit vectors, the cosine similarity is equal to the vector dot product.

[0061]

[0062] in, For pixels The similarity value ranges from [0,1]. A higher value indicates that the pixel is more similar to the rubber plantation sample in terms of features. is the feature dimension of the satellite embedded image (SE), with a maximum value of 64; j is the dimension index; For the first The first pixel Dimensional embedding feature values.

[0063] 3. Similarity layer generation: Combine the similarity values ​​of all pixels to generate a single-band similarity image, which has the same spatial resolution and georeferenced information as the original embedded image.

[0064] Step S4: Low-confidence cell removal based on dynamic threshold

[0065] 1. Obtain the initial classification result image to be corrected and convert it into a binary image, where the target land cover category pixel is assigned a value of 1 and other category pixels are assigned a value of 0.

[0066] 2. Construct a spatial neighborhood analysis kernel. Define the kernel radius r (typically r=1), generating a square convolution kernel of size (2r+1)×(2r+1). Set the weight at the kernel center position to 0 (excluding the center pixel itself), and set the weights at other positions to 1. The maximum possible number of neighboring pixels within the neighborhood ( )for:

[0067]

[0068] 3. Perform neighborhood convolution on the binary classification image, count the number of pixels in the neighborhood of each pixel that are classified as the target feature, and generate a neighborhood counting image. The maximum number of neighbors within the neighborhood window.

[0069] 4. Determine the similarity threshold. The threshold is dynamically calculated using the percentile method:

[0070] 4.1 Construct an effective pixel mask: Select pixels that simultaneously meet the two conditions of "already classified as the target ground feature" and "the number of neighbors is greater than or equal to the minimum neighborhood threshold (removed)".

[0071] 4.2 Extract the similarity values ​​corresponding to the effective pixels, and calculate the p-th percentile as the removal threshold θ (typical values ​​p=1, 2, or 5):

[0072]

[0073] In the formula, The dataset consists of the similarity values ​​of these valid pixels. Percentiles specified by the user.

[0074] This dynamic thresholding method allows the removal intensity to adapt to the data distribution, avoiding the problem of inconsistent performance of a fixed threshold in different regions.

[0075] 5. Perform pixel removal determination. For each pixel already classified as a target feature, determine whether to retain it based on the following logic:

[0076] • If the number of neighbors is greater than or equal to the minimum neighborhood threshold and the similarity is greater than or equal to the dynamic threshold θ, then the cell is retained.

[0077] • If the number of neighbors is less than the minimum neighborhood threshold (to remove), the cell is retained to prevent valid isolated patches from being incorrectly removed;

[0078] • In other cases, remove the pixel;

[0079] 6. Generate a refined classification image after removal, and record the spatial distribution of the removed pixels for subsequent analysis.

[0080] Step S5: Gap filling based on neighborhood density and similarity

[0081] 1. Identify gap cells. Compare the refined classified image with the original classified image and mark the locations of all cells not classified as target features as candidate filling areas.

[0082] 2. Construct the padding analysis kernel. Define the padding kernel radius R (typically R=2, larger than the kernel radius of the removal stage to obtain a larger spatial context), and generate the convolution kernel using the same method as in step 4. The padding stage uses a larger neighborhood range to ensure that the padding decision is based on sufficient spatial evidence.

[0083] 3. Calculate neighborhood statistics:

[0084] 3.1 Neighborhood Counting: Counting the number of pixels in the neighborhood of each gap cell that have been classified as the target land cover. .

[0085] 3.2 Neighborhood Coverage: Calculate the proportion of classified neighborhoods to the maximum possible number of neighborhoods.

[0086]

[0087] In the formula, This represents the number of pixels in the neighborhood that have already been classified as the target category. It is the maximum possible number of neighbors within this neighborhood window.

[0088] 3.3 Average Neighborhood Similarity: Calculate the average embedding similarity between the gap cell and its classified neighbors. This calculation requires extracting the embedding vector at each neighbor location, calculating the dot product with the center cell, and averaging over all classified neighbors.

[0089]

[0090] in, This represents the average similarity value of the classified neighborhoods within the gap cell's neighborhood. For the central pixel and the first Embedding similarity of neighborhood locations, For the first The classification status (0 or 1) of each neighborhood location.

[0091] 4. Perform filling determination. A gap cell must meet all of the following conditions to be filled as the target feature:

[0092] • Condition 1: Number of neighbors ≥ Minimum neighborhood threshold (filling) (typical value 8)

[0093] • Condition 2: Neighborhood coverage rate ≥ minimum coverage rate threshold (typical value 0.3)

[0094] • Condition 3: Average neighborhood similarity ≥ minimum similarity threshold (typical value 0.85)

[0095] 5. Optionally, additional constraint masks (such as slope masks, elevation masks, water body masks, etc.) can be superimposed to further limit the filling range and prevent accidental filling in obviously unsuitable areas.

[0096] 6. Assign the target land cover category to the pixels that meet the filling conditions, merge them with the refined classification image, and generate the final corrected classification result image.

[0097] Step S6: Output and Evaluation

[0098] 1. Output the corrected classification result image, maintaining the same spatial reference and data format as the input image.

[0099] 2. Optional output of intermediate diagnostic images, including: similarity images, neighborhood count images, distribution of removed pixels, distribution of filled pixels, etc., for correction effect analysis and parameter tuning.

[0100] 3. Calculate classification accuracy indices (overall accuracy, Kappa coefficient, etc.) before and after correction using independent validation samples to quantitatively evaluate the correction effect.

[0101] The key parameters are explained below:

[0102] Example

[0103] Parameter settings: Percentile: 1; Denoising kernel radius: 1; Filling kernel radius: 2; Minimum number of neighbors: 5; Minimum coverage: 0.60; Minimum similarity: 0.90. The results are as follows. Figure 1 As shown: a) is a high-resolution Google base map of the study area; b) is the original classification and distribution results of rubber forests in the study area, with red areas representing rubber forests; c) is the final classification and distribution results of rubber forests after correction with embedded satellite data; d) is a magnified view of a local sample area; e) is the original classification and distribution results of rubber forests in the corresponding local sample area; f) is the final classification and distribution results of rubber forests in the corresponding local sample area after correction with embedded satellite data.

[0104] Taking the 2024 classification results of rubber plantations on Hainan Island as an example, after independent validation of the sample set, the overall accuracy of the rubber plantation classification results optimized by the method of this invention reached 92.20%, with a Kappa coefficient of 0.8428; while the overall accuracy of the unoptimized initial classification results was 90.89%, with a Kappa coefficient of 0.8163. The comparison shows that this method improves the overall accuracy by 1.31% and the Kappa coefficient by 0.0265, effectively solving problems such as fragmented patches, internal voids, and blurred edges in the initial classification results. It significantly improves the spatial consistency and feature matching degree of rubber plantation classification, demonstrating a significant optimization effect and providing high-confidence data support for rubber plantation resource inventory, dynamic monitoring, and refined management.

[0105] Comparative Example

[0106] Taking rubber plantation classification as an example, the method of this invention is used to correct the initial classification results. The accuracy comparison under different removal intensities is shown in the table below:

[0107]

[0108] Experimental results show that the method of the present invention can effectively improve classification accuracy. With the configuration of removing 1% of low confidence pixels, the overall accuracy is improved by 1.31 percentage points and the Kappa coefficient is improved by 0.0265.

[0109] Although embodiments of the invention have been shown and described, it will be understood by those skilled in the art that various changes, modifications, substitutions and alterations can be made to these embodiments without departing from the principles and spirit of the invention, the scope of which is defined by the appended claims and their equivalents.

Claims

1. An optimization method for rubber plantation classification results based on satellite embedded (SE) features, characterized in that, Includes the following steps: S1. Satellite Embedded Image Data Acquisition and Processing: Acquire satellite embedded images covering the target area, wherein each pixel of the satellite embedded image contains an N-dimensional embedding vector; S2. Construction of embedded feature template for rubber plantation: Based on the known rubber plantation reference sample area, the embedding vector of all pixels in the sample area is extracted, the average value is calculated and normalized to obtain the feature template vector of the rubber plantation. S3. Calculate the similarity of pixels in the entire region: Traverse each pixel in the target region, calculate the cosine similarity between its embedding vector and the feature template vector, and generate a similarity layer for the entire region. S4. Low-confidence pixel removal based on dynamic threshold: Based on the initial rubber plantation classification results, combined with the spatial neighborhood density of the pixels and their similarity values ​​in the similarity layer, a dynamically calculated percentile threshold is used to identify and remove low-confidence pixels in the classification results to obtain a refined classification image. S5. Gap filling based on neighborhood density and spectral similarity: Identify gap pixels in the refined classification image, fill gap pixels that simultaneously meet the conditions of the number of classified pixels in the neighborhood, neighborhood coverage, and neighborhood average spectral similarity threshold, and generate the final optimized rubber plantation classification result.

2. The method for optimizing rubber plantation classification results based on satellite embedded (SE) features according to claim 1, characterized in that, The specific process of constructing the rubber plantation embedded feature template in step S2 includes: S21. Generate a binary sample mask image based on the rubber plantation reference sample area; S22. Overlay the satellite embedded image with the sample mask image and extract the embedding vector of all pixels in the sample area; S23. Calculate the average value of all pixel embedding vectors in the sample region to obtain the initial feature template vector; S24. Perform L2 normalization on the initial feature template vector to obtain the feature template vector in unit vector form.

3. The method for optimizing rubber plantation classification results based on satellite embedded (SE) features according to claim 1, characterized in that, In step S3, the cosine similarity algorithm is used to calculate the similarity value between the embedding vector of each pixel and the feature template vector. The calculation formula is as follows: in, For pixels The similarity value ranges from [0,1]. A higher value indicates that the pixel is more similar to the rubber plantation sample in terms of features. is the feature dimension of the satellite embedded image (SE), with a maximum value of 64; j is the dimension index; For the first The first pixel Dimensional embedding feature values.

4. The method for optimizing rubber plantation classification results based on satellite embedded (SE) features according to claim 1, characterized in that, The specific process of low-confidence pixel removal based on dynamic threshold in step S4 includes: S41. Convert the initial rubber plantation classification results into binary images; S42. Perform neighborhood analysis on the binary image and count the number of pixels that have been classified as rubber plantations within the predetermined neighborhood of each pixel; S43. Construct a set of effective pixels, wherein the effective pixels are pixels that have been classified as rubber plantations and whose number of neighbors is greater than or equal to the minimum number of neighbors (removed); S44. Calculate the p-th percentile of the similarity value corresponding to the effective pixel set in the similarity layer, and use it as the dynamic removal threshold θ; S45. For each pixel classified as a rubber plantation, if its number of neighbors is less than the minimum number of neighbors (removed), it is retained; if its number of neighbors is greater than or equal to the minimum number of neighbors (removed) and its similarity value is greater than or equal to the dynamic removal threshold θ, it is retained; otherwise, it is removed from the classification results.

5. The method for optimizing rubber plantation classification results based on satellite embedded (SE) features according to claim 1, characterized in that, The specific process of filling the gaps in step S5 includes: S51. Compare the refined classification image with the initial classification result and identify pixels that were not classified as rubber plantations as candidate gap-filling pixels. S52. For each of the candidate filling gap pixels, the following judgment is made: Calculate the number of pixels within its predetermined filling neighborhood that have been classified as rubber plantations. ; Calculate the neighborhood coverage rate, i.e. The ratio of the number of possible neighbors to the number of neighbors within the filled neighborhood; Calculate the average embedding similarity value between the gap cell and all its neighboring cells that have been classified as rubber plantations. ; S53. If the candidate filling gap pixel simultaneously meets the following conditions, it shall be classified as a rubber plantation; (a) ≥ Minimum neighborhood number (filled); (b) The neighborhood coverage rate is ≥ the coverage threshold; (c) ≥ Similarity threshold.

6. The method for optimizing rubber plantation classification results based on satellite embedded (SE) features according to claim 4, characterized in that, The percentile p in step S44 ranges from 1 to 5.

7. The method for optimizing rubber plantation classification results based on satellite embedded (SE) features according to claim 5, characterized in that, In step S52, the range of the neighborhood to be filled is greater than the range of the neighborhood analysis in step S42 during the removal phase.

8. The method for optimizing rubber plantation classification results based on satellite embedded (SE) features according to claim 5, characterized in that, Following step S53, the method further includes the step of: limiting the filling result using an additional geographic constraint mask, the geographic constraint mask including at least one of a slope mask, an elevation mask, or a water body mask.