A weakly supervised detection method for bulldozing areas based on high-resolution remote sensing images

Abstract

The application provides a kind of weak supervision detection method of bulldozing area based on high-resolution remote sensing image, the method comprises: the original remote sensing image is preprocessed and feature extraction, obtain tile set and corresponding multi-scale feature map;Seed point mining is carried out to multi-scale feature map, and the seed point set of bulldozing area is obtained;Based on the seed point set of bulldozing area, the seed point corresponding to the multi-scale feature map is determined, and the complete region activation map is generated by expanding from the seed point;Based on the original remote sensing image and multi-scale feature map, the complete region activation map is processed by multi-feature boundary perception algorithm, and the pseudo label corresponding to the boundary position is obtained;Based on tile set and pseudo label, the new remote sensing image is processed by semantic segmentation network, and the detection result corresponding to the new remote sensing image of bulldozing area is obtained.The application realizes land use survey, illegal land use law enforcement inspection and land change dynamic monitoring by only relying on pixel level bulldozing area accurate detection of image level label.

Description

Technical Field

[0001] This invention relates to the field of remote sensing image interpretation technology, and in particular to a method for weakly supervised detection of bulldozing areas based on high-resolution remote sensing images. Background Technology

[0002] With the rapid development of Earth observation technology, high-resolution remote sensing imagery provides an unprecedented data foundation for monitoring surface changes, especially for detecting drastic surface disturbances caused by human activities (such as bulldozing areas). Bulldozing areas refer to bare surface areas formed after the removal of existing vegetation and cover due to engineering construction, land reclamation, infrastructure construction, or illegal occupation. As a direct indicator of engineering construction and land leveling activities, rapid and accurate automated detection of bulldozing areas is of great significance for urban planning, land monitoring, ecological environment assessment, and disaster emergency response.

[0003] Currently, automated bulldozing detection technology based on high-resolution remote sensing imagery faces a long-standing dilemma: balancing application cost, method adaptability, and feature representation. High-precision fully supervised deep learning methods heavily rely on massive amounts of pixel-level labeled data, a process that is extremely time-consuming and labor-intensive, leading to high model training costs and hindering the support of large-scale, dynamic real-world monitoring needs. Weakly supervised learning methods introduced to reduce labeling dependence, such as category activation map (CAM) techniques, encounter significant challenges in the specific scenario of bulldozing areas: bulldozing areas typically appear as small-scale targets with a very low proportion in images, resulting in weak and fragmented initial activation responses generated by traditional methods, making it difficult to fully locate the target and leading to severe missed detections.

[0004] Furthermore, the inherent complexity and high intra-class variability of bulldozing areas further exacerbate the identification difficulty. Specifically, the spectral characteristics of bulldozing areas are easily confused with those of bare soil, sandy land, and other features. Moreover, the morphology, scale, and texture vary greatly due to differences in construction stage and scale, resulting in insufficient generalization ability of models relying on single features or fixed scales. Therefore, existing technologies often fall into a series of bottlenecks when dealing with bulldozing area detection tasks, including prohibitive annotation costs, failure of weakly supervised methods, and difficulty in robustly representing complex features. An innovative method that can simultaneously solve these problems is urgently needed. Summary of the Invention

[0005] The purpose of this invention is to provide a weakly supervised detection method for bulldozing areas based on high-resolution remote sensing images, which solves the shortcomings of existing methods in small target detection, feature representation, and pseudo-label quality, and achieves high-precision pixel-level bulldozing area segmentation under the condition of only image-level labels, while reducing labeling costs and maintaining detection accuracy close to that of fully supervised methods.

[0006] To achieve the above objectives, this invention provides a method for weakly supervised detection of bulldozing areas based on high-resolution remote sensing images, the method comprising: S11. Preprocess and extract features from the original remote sensing images to obtain a set of patches and corresponding multi-scale feature maps; S12. Perform seed point mining on the multi-scale feature map to obtain the seed point set of the bulldozing area; S13. Based on the seed point set of the bulldozing area, determine the seed points corresponding to the multi-scale feature map, and expand from the seed points to generate a complete region activation map. S14. Based on the original remote sensing image and multi-scale feature map, the activation map of the complete region is processed by the multi-feature boundary perception algorithm to obtain the pseudo label corresponding to the boundary position. S15. Based on the tile set and pseudo-labels, the new remote sensing image is processed by a semantic segmentation network to obtain the bulldozing area detection results corresponding to the new remote sensing image.

[0007] Furthermore, a set of map tiles and corresponding multi-scale feature maps are generated, specifically including: S21. Preprocess the original remote sensing image to obtain a set of tiles, perform binary classification and labeling on the set of tiles to form a tile dataset; S22. Construct a multi-scale feature pyramid for each tile in the tile dataset and output the target image; S23. Extract features from the target image and output a multi-scale feature map.

[0008] Furthermore, seed point mining is performed on the multi-scale feature map, specifically including: S31. Based on the historical background category diagram, construct a background category prototype library; S32. Extract the feature vector of each spatial location in the multi-scale feature map, calculate the Euclidean distance between the feature vector and all backgrounds in the background category prototype library, and take the minimum distance as the background belonging degree of the corresponding location in the multi-scale feature map. S33. Convert the background attribution degree into a contrast significance value, and use a Gaussian function for nonlinear mapping to obtain a contrast significance map; S34. The contrast saliency maps are fused to obtain a fused saliency map. Thresholding and local maximum detection are performed on the fused saliency map to extract the seed point set of the bulldozing area.

[0009] Furthermore, generating a complete region activation map specifically includes: S41. Flatten the multi-scale feature map into several feature vectors, calculate the semantic similarity between any two feature vectors, and construct the semantic affinity matrix between pixels; S42. Normalize the semantic affinity matrix to obtain the transition probability matrix, encode the seed point set of the bulldozing area into the initial activation vector, and generate the activation region. S43. Based on the transition probability matrix, the activation region is expanded by iterative matrix multiplication, characterized in that a complete activation map of the region is obtained.

[0010] Furthermore, a spatial distance penalty is added to the calculation of the semantic similarity between any two feature vectors.

[0011] Furthermore, generating pseudo-labels corresponding to the boundary positions specifically includes: S51. Based on the preprocessed original remote sensing image, calculate the gradient magnitude of each band and fuse them to obtain the spectral gradient boundary map. S52. Extract multi-directional and multi-scale texture features from the original remote sensing image to obtain a texture boundary map; S53. Extract the target boundary information from the multi-scale feature map to obtain the depth feature boundary map; S54. Fuse the spectral gradient boundary map, texture boundary map, and depth feature boundary map to generate a fused boundary map; S55. Optimize the activation map of the complete region by fusing the boundary map to obtain the pseudo-labels corresponding to the boundary positions.

[0012] Furthermore, the new remote sensing images are processed using a semantic segmentation network, specifically including: S61. Preprocess and extract features from the new remote sensing image to obtain a set of map patches. Process the set of map patches through a trained semantic segmentation network to output a probability map of the bulldozing area. S62. After converting the bulldozing area probability map into a binary mask, the pixels in the binary mask that are greater than a preset threshold are marked as bulldozing areas to obtain the bulldozing area pixel map. S63. After performing opening and closing operations on the pixel image of the bulldozing area, remove the patches with an area smaller than the preset threshold. S64. Convert the binary mask in the remaining patches into vector polygons, calculate the geometric properties of each patch, and obtain the bulldozing area map; S65. Assign operational attributes to each bulldozing area patch in the bulldozing area map to obtain the bulldozing area detection results corresponding to the new remote sensing image.

[0013] Furthermore, each bulldozing area patch in the bulldozing area map is assigned business attributes, specifically including: S71. Obtain at least one first set of geographic features, wherein the first set of geographic features contains geospatial layer data with legal management attributes; S72. Perform spatial overlay analysis on each bulldozed area map patch and at least one set of first geographic elements, and assign at least one legal management attribute information from the set of first geographic elements to each bulldozed area map patch based on the analysis results; S73. Based on preset rules, conduct compliance judgment on bulldozing area patches that have been assigned legal management attribute information, and generate bulldozing area detection results that include legal management attribute information and compliance judgment results.

[0014] Compared with the prior art, the beneficial effects of the present invention are: This invention provides a weakly supervised detection method for bulldozing areas based on high-resolution remote sensing images. Firstly, traditional CAM methods directly utilize the gradient information of classification networks to generate activation maps, resulting in weak responses to small targets. The contrast saliency seed point mining mechanism proposed in this invention amplifies the discriminative power of small targets at the feature level by constructing a contrast feature space with the background category. Experiments show that the recall rate for bulldozing areas accounting for less than 3% is increased from 45% in traditional CAM to 78%. Secondly, existing weakly supervised methods are mainly designed for salient targets in natural images, without considering the characteristics of remote sensing images. This invention is specifically designed for bulldozing area detection scenarios, integrating spectral, texture, and depth features for boundary optimization, significantly enhancing the ability to distinguish bulldozing areas from easily confused features such as bare soil and sand, improving segmentation IoU by 8-12 percentage points.

[0015] Finally, while fully supervised methods require pixel-level annotation, this invention only requires image-level labels, reducing annotation costs by over 95%. It achieves 85%-90% of the performance of fully supervised methods while significantly lowering the barrier to practical application. For data from new regions and time periods, annotation and model adaptation can be completed quickly, supporting large-scale, routine monitoring operations. This invention proposes a small-target perception progressively extended network to achieve pixel-level accurate detection of bulldozing areas relying solely on image-level labels, providing a highly efficient and intelligent technical means for natural resources departments to conduct land use surveys, illegal land use enforcement inspections, and dynamic monitoring of land changes. Attached Figure Description

[0016] To more clearly illustrate the technical solutions in the embodiments of this application or the prior art, the drawings used in the description of the embodiments or the prior art will be briefly introduced below. Obviously, the drawings described below are only embodiments of this application. For those skilled in the art, other drawings can be obtained based on the provided drawings without creative effort.

[0017] Figure 1 This is a schematic diagram of a method for weakly supervised detection of bulldozing areas based on high-resolution remote sensing images, provided in an embodiment of the present invention. Figure 2This is a schematic diagram of the overall framework of a weakly supervised detection method for bulldozing areas based on high-resolution remote sensing images, provided in an embodiment of the present invention. Figure 3 A schematic diagram illustrating the principle of comparative saliency seed point mining and progressive expansion provided in an embodiment of the present invention; Figure 4 This is a schematic diagram illustrating the principle of multi-feature boundary sensing fusion provided in an embodiment of the present invention. Detailed Implementation

[0018] The present invention will now be described in further detail with reference to the accompanying drawings and embodiments. It should be understood that the specific embodiments described herein are merely illustrative of the invention and not intended to limit it. Furthermore, it should be noted that, for ease of description, the accompanying drawings show only the parts relevant to the present invention, and not all of the structures.

[0019] Reference Figure 1 and Figure 2 This embodiment provides a weakly supervised detection method for bulldozing areas based on high-resolution remote sensing images. The method includes: S11. Preprocess and extract features from the original remote sensing images to obtain a set of patches and corresponding multi-scale feature maps.

[0020] S12. Perform seed point mining on the multi-scale feature map to obtain the seed point set of the bulldozing area; S13. Based on the seed point set of the bulldozing area, determine the seed points corresponding to the multi-scale feature map, and expand from the seed points to generate a complete region activation map.

[0021] S14. Based on the original remote sensing image and multi-scale feature map, the activation map of the complete region is processed by the multi-feature boundary perception algorithm to obtain the pseudo-labels corresponding to the boundary positions.

[0022] S15. Based on the tile set and pseudo-labels, the new remote sensing image is processed by a semantic segmentation network to obtain the bulldozing area detection results corresponding to the new remote sensing image.

[0023] In this embodiment, the first layer involves seed point mining. Addressing the weak CAM response for small targets, a Contrastive Saliency Seed Mining (CSSM) mechanism is proposed. By constructing a contrast feature space between the bulldozing area and background features, the discriminative power of small targets is amplified at the feature level, reliably locating the initial seed points of the bulldozing area. The second layer involves region expansion, proposing a Semantic Consistency Guided Progressive Expansion (SCPE) module. Starting from the seed points, based on semantic similarity and spatial continuity constraints between pixels, the activation region is gradually expanded to the complete bulldozing area, solving the problem that CAM can only locate discriminative regions and cannot obtain the complete outline. The third layer involves boundary optimization, designing a Multi-feature Boundary Awareness (MBA) module that fuses spectral gradient, texture edge, and depth feature boundary information to refine the boundary position of the bulldozing area and generate high-quality pseudo-labels for training the final segmentation network.

[0024] As a preferred embodiment, generating a set of tiles and corresponding multi-scale feature maps specifically includes: S21. Preprocess the original remote sensing image to obtain a set of tiles, perform binary classification and labeling on the set of tiles to form a tile dataset.

[0025] S22. Construct a multi-scale feature pyramid for each tile in the tile dataset and output the target image.

[0026] S23. Extract features from the target image and output a multi-scale feature map.

[0027] In this embodiment, the original remote sensing image is preprocessed and its features are extracted to provide multi-scale, multi-level feature representations for subsequent seed point mining and region expansion. The input data is high-resolution multispectral remote sensing imagery (such as data from domestic high-resolution satellites like Gaofen-2 (1-meter resolution), Gaofen-1 (2-meter resolution), and Beijing-3). The preprocessing workflow includes: (1) Radiometric calibration: Convert the raw digital quantization values ​​(DN values) recorded by the sensor into radiance values ​​or reflectance values ​​with clear physical meaning, eliminate system differences between different sensors, and make the image data comparable; (2) Atmospheric correction: The FLAASH or 6S atmospheric radiative transfer model is used to eliminate the influence of aerosol scattering and water vapor absorption in the atmosphere on the surface reflection signal, obtain the true surface reflectance, and improve the consistency of images at different time phases. (3) Orthorectification: Using digital elevation model (DEM) and satellite orbit parameters, geometric distortions caused by terrain undulations and sensor attitude changes are eliminated, so that the image has the correct geographic coordinates; (4) Image Tile Segmentation: The large image processed as described above is segmented into fixed-size tiles (512×512 pixels in this invention). Adjacent tiles maintain a 10% overlap (approximately 51 pixels) to avoid truncating bulldozing targets located at tile boundaries and affecting detection performance. The segmented tile set is denoted as {I_1, I_2, ..., I_N}, where N is the total number of tiles.

[0028] Each segmented image patch is binary-classified and labeled with a label y∈{0,1}. Here, y=1 indicates that the patch contains a bulldozing area (positive sample), and y=0 indicates that the patch does not contain a bulldozing area (negative sample). The labeling does not require specifying the exact location, shape, or extent of the bulldozing area; it only needs to determine whether the patch as a whole contains a bulldozing area. This image-level labeling method significantly reduces the labeling difficulty. A skilled labeler can complete the labeling of 3000-5000 images per day, which is 50-100 times more efficient than pixel-level labeling. After labeling, a labeled image patch dataset {(I_n, y_n)} is formed, where n=1,2,...,N.

[0029] Because the scale of targets in bulldozing areas varies greatly (from tens of square meters to several hectares), features at a single scale are insufficient to effectively represent all targets. Therefore, a multi-scale feature pyramid is constructed for the patch dataset. This increases the relative proportion of small targets at smaller scales, resulting in stronger feature responses, while preserving spatial details at larger scales facilitates subsequent boundary refinement. Let the input patch be I∈R^(H×W×C), where H represents the height of the patch (number of pixels), W represents the width of the patch (number of pixels), and C represents the number of spectral channels (e.g., C=3 for RGB three-channel images, C=4 or more for multispectral images). An image pyramid containing L scales is constructed, as shown in the following formula: I_l = Resize(I, s_l),l = 1, 2, ..., L Where I_l represents the scaled image at the l-th scale, i.e., the target image; Resize(·) represents the bilinear interpolation scaling operation; s_l represents the scaling factor at the l-th scale, with a value range of (0,1). This invention uses L=4 scales, with scaling factors of s={1.0, 0.75, 0.5, 0.25}, i.e., the original size, 3 / 4 size, 1 / 2 size, and 1 / 4 size. By constructing a multi-scale pyramid, the relative area proportion of a small bulldozing target in a small-scale image (e.g., s=0.25) is increased by more than 4 times, thereby obtaining a stronger response signal in subsequent feature extraction. For example, a small bulldozing area that occupies only 1% of the area in the original image will increase its relative proportion to about 16% in the 1 / 4 scale image, significantly improving the detectability of small targets.

[0030] For the target image, a pre-trained visual backbone network is used to extract deep features. This invention supports ResNet-50 or Vision Transformer (ViT) as the backbone network. Taking ResNet-50 as an example, feature maps output from the three convolutional stages conv3, conv4, and conv5 are extracted. The spatial resolutions of these three stages correspond to 1 / 8, 1 / 16, and 1 / 32 of the original image, respectively, forming a multi-level feature representation from shallow to deep layers and from details to semantics. F_l = {F_l^3, F_l^4, F_l^5},F_l^k ∈ R^(H_k×W_k×D_k) Where F_l represents the multi-level feature set of the image at scale l, i.e., the multi-scale feature map; F_l^k represents the feature map of the k-th layer at scale l (k=3,4,5 corresponding to conv3, conv4, conv5); H_k and W_k represent the height and width (in pixels) of the k-th layer feature map, respectively; D_k represents the channel dimension of the k-th layer feature map (512 dimensions for conv3, 1024 dimensions for conv4, and 2048 dimensions for conv5). The multi-scale image pyramid combined with multi-level feature extraction produces a total of L×3=12 sets of feature maps. Shallow features (conv3) contain more texture and edge details, which is beneficial for boundary localization; deep features (conv5) contain stronger semantic information, which is beneficial for object recognition. These multi-scale, multi-level features will be passed to the subsequent Contrast Saliency Seed Mining (CSSM) module for seed point mining.

[0031] As a preferred embodiment, seed point mining is performed on the multi-scale feature map, specifically including: S31. Based on the historical background category diagram, construct a background category prototype library.

[0032] S32. Extract the feature vector of each spatial location in the multi-scale feature map, calculate the Euclidean distance between the feature vector and all backgrounds in the background category prototype library, and take the minimum distance as the background belonging degree of the corresponding location in the multi-scale feature map.

[0033] S33. Convert the background attribution degree into a contrast significance value, and use a Gaussian function for nonlinear mapping to obtain a contrast significance map.

[0034] S34. The contrast saliency maps are fused to obtain a fused saliency map. Thresholding and local maximum detection are performed on the fused saliency map to extract the seed point set of the bulldozing area.

[0035] In this embodiment, refer to Figure 3 The upper half of the diagram (CSSM principle) uses an abstract feature space diagram to illustrate the response distribution of traditional CAM in the feature space. The left side shows the response of small targets (small dots) in the feature space, with weak responses. The right side shows the contrast saliency method, which, by introducing background category prototypes (different shaped markers), makes the distance between small targets and the background explicit, thus enhancing the saliency response. The lower half (SCPE principle) uses a grid diagram to illustrate the progressive expansion process. The first frame shows the initial position of the seed point (black square) in the grid; subsequent frames show the process of the activation region (grayscale representing activation intensity) gradually expanding outward; the last frame shows the complete target region after expansion, with arrows indicating the iteration direction.

[0036] Based on multi-scale feature maps, this invention identifies and locates high-confidence initial seed points in bulldozing areas by constructing a contrastive feature space with background features. The contrastive saliency seed point mining module is the first core innovation of this invention, aiming to solve the problem of weak CAM response for small targets. The core idea is: instead of directly generating class activation maps for bulldozing areas (traditional methods have weak responses to small targets), it constructs a contrastive feature space between bulldozing areas and various background features, calculating the feature distance between each pixel location and known background categories. The larger the distance, the less likely the location is to be background, i.e., the more likely it is to be a bulldozing area. Through this contrast mechanism, the discriminative power of small targets is significantly amplified at the feature level.

[0037] First, a category prototype library covering the main background land cover types in the study area is constructed. The category prototype library is denoted as B = {b_1, b_2, ..., b_M}, where M is the total number of background categories. Typical background categories include: vegetation (farmland, woodland, grassland, orchards); artificial surfaces (buildings, roads, paved surfaces); water bodies (rivers, lakes, ponds); and bare land (natural bare soil, sandy land, Gobi desert) (categories easily confused with bulldozed areas require careful differentiation). Representative samples for each background category are collected, utilizing existing land use survey data from the natural resources department (such as the Third National Land Survey data), extracting typical map patches for each category as samples. 50-100 sample patches are collected for each category. For each sample patch, a backbone network (ResNet-50) is used to extract depth features. For pixel regions belonging to that category within the sample patch, the corresponding feature vector is extracted. The mean of all sample feature vectors for the m-th background category is calculated to obtain the prototype feature vector for that category. The calculation formula is as follows: p_m = (1 / N_m) × Σ_{i=1}^{N_m} f(x_i^m) Where p_m represents the prototype feature vector of the m-th background feature, which is a D-dimensional vector (D is the feature dimension, such as 2048); N_m represents the total number of sample pixels of the m-th background; x_i^m represents the i-th sample pixel of the m-th feature; f(·) represents the feature extraction function of the backbone network, which maps the input pixels to a D-dimensional feature vector; Σ represents the summation of all sample feature vectors. By averaging a large number of samples, prototype vectors that can represent the typical characteristics of the category are obtained. The prototype vectors eliminate the random fluctuations of individual samples and retain the common features of the category. The final constructed category prototype library contains M prototype vectors {p_1, p_2, ..., p_M} for subsequent comparison significance calculation.

[0038] Using a class prototype library, the contrast saliency is calculated for each spatial location in the multi-scale feature map. The feature vector f_{i,j} for each spatial location (i,j) in the input multi-scale feature map is extracted. The Euclidean distance between the feature vector and all M background class prototypes is calculated, and the minimum distance is taken as the background attribution. d_bg(i,j) = min_{m∈{1,2,...,M}} ||f_{i,j} - p_m||_2 Where d_bg(i,j) represents the background affiliation score at position (i,j), which is a non-negative real number; min represents the minimum value operation; m∈{1,2,...,M} represents traversing all M background categories; f_{i,j} represents the feature vector at position (i,j); p_m represents the prototype vector of the m-th background category; ||·||_2 represents the L2 norm (Euclidean distance). The background affiliation score d_bg reflects the feature difference between the current position and the nearest background category. If the d_bg value is small, it means that the features at this position are similar to a known background category, and it is very likely to belong to that background category; if the d_bg value is large, it means that the features at this position are not similar to any known background categories, and it may be a new category not covered by the background library, i.e., a potential bulldozer zone.

[0039] The background attribution is converted into a contrast significance value, and a Gaussian kernel function is used for nonlinear mapping, as shown in the following formula: S_contrast(i,j) = 1 - exp(-d_bg(i,j)² / (2σ²)) Where S_contrast(i,j) represents the contrast significance value of position (i,j), with a value range of [0,1); exp(·) represents the natural exponential function; d_bg(i,j)² represents the square of the background attribution; σ is the bandwidth parameter of the Gaussian kernel, controlling the sensitivity of the saliency response. In this invention, σ is taken as the median of the background attribution for all positions. When d_bg is close to 0 (similar to the background), S_contrast is close to 0, indicating that the position is not significant; when d_bg is large (largely different from the background), S_contrast is close to 1, indicating that the position is highly significant. By calculating the saliency of all positions in the multi-scale feature map, a contrast significance map S_contrast∈R^(H×W) with the same spatial resolution as the original image is obtained. The highlighted areas in the map are the potential bulldozing areas. Thus, even if the bulldozing area is small, as long as its features differ from all background categories, it will produce a significant response, thereby solving the problem of weak response of traditional CAM to small targets.

[0040] For each scale of the multi-scale feature map, the saliency is calculated, resulting in L saliency maps at different scales. Since the spatial resolution of the images at each scale is different, they need to be uniformly upsampled to the original image size before fusion. S_fused = Σ_{l=1}^L w_l × Upsample(S_contrast^l) Wherein, S_fused represents the fused saliency map, i.e., the fused saliency map, with a size of H×W; Σ represents the weighted sum of L scales; w_l represents the fusion weight of the l-th scale, satisfying Σw_l=1, and this invention adopts equal weights w_l=1 / L; Upsample(·) represents the bilinear interpolation upsampling operation, which uniformly adjusts the saliency maps of each scale to the original image size H×W; S_contrast^l represents the contrast saliency map of the l-th scale. Multi-scale fusion integrates target response information at different scales: small scales (e.g., s=0.25) are sensitive to small targets, while large scales (e.g., s=1.0) retain more spatial details, and the fused saliency map can simultaneously detect bulldozing targets of different sizes.

[0041] Thresholding and local maxima detection are performed on the fusion saliency map to extract the seed point set: Seeds = {(i,j) | S_fused(i,j)>τ_high and (i,j) is a local maximum point} Here, Seeds represents the extracted set of seed points, i.e., the bulldozing zone seed point set, where each element is a two-dimensional coordinate (i,j); τ_high is the high confidence threshold, which in this invention is taken as the 90th percentile of the numerical distribution of the fused saliency map, meaning that only positions ranking in the top 10% of saliency are likely to become seed points; local maxima points refer to the positions with the largest saliency values ​​within a 3×3 neighborhood, used to exclude redundant responses. Through dual screening using the high threshold and local maxima, the extracted seed points are ensured to have high confidence and accurate location. Each seed point represents a location that is highly likely to be a bulldozing zone and will serve as the starting point for the next stage of semantic consistency progressive expansion module (SCPE) region expansion. The seed point set Seeds is the final output of this stage.

[0042] As a preferred embodiment, generating a complete region activation map specifically includes: S41. Flatten the multi-scale feature map into several feature vectors, calculate the semantic similarity between any two feature vectors, and construct the semantic affinity matrix between pixels.

[0043] S42. Normalize the semantic affinity matrix to obtain the transition probability matrix, encode the seed point set of the bulldozing area into the initial activation vector, and generate the activation region.

[0044] S43. Based on the transition probability matrix, the activation region is expanded by iterative matrix multiplication, characterized in that a complete activation map of the region is obtained.

[0045] As another preferred embodiment, a spatial distance penalty is added to the calculation of the semantic similarity between any two feature vectors.

[0046] In this embodiment, starting with seed points in the seed point set, region expansion is performed on the multi-scale feature map, expanding isolated seed points into activated regions covering the entire bulldozing area. The semantically consistent progressive expansion module is the second core innovation of this invention, solving the problem of expanding from discrete seed points to the complete target region. That is, starting with seed points, based on the semantic similarity between pixels (the degree of similarity between feature vectors), the activated region is iteratively expanded outwards, activating neighboring pixels semantically similar to the seed points, until the expansion automatically stops at the target boundary (where a semantic abrupt change occurs). This process simulates the classic region growing algorithm, but is performed in a deep feature space, resulting in stronger semantic consistency.

[0047] Based on the multi-scale feature map (using features from the conv4 layer, with a resolution of 1 / 16 of the original image), a semantic affinity matrix between pixels is constructed. The multi-scale feature map F∈R^(H'×W'×D) is flattened into N=H'×W' feature vectors, and the semantic similarity between any two positions is calculated: A(i,j) = exp(f_i^T × f_j / (||f_i|| × ||f_j|| × τ)) Where A(i,j) represents the semantic affinity (similarity) between positions i and j, with a value range of (0,+∞); f_i and f_j represent the feature vectors (D-dimensional) of positions i and j, respectively; f_i^T represents the transpose of f_i; f_i^T×f_j represents the inner product of the two vectors, reflecting the consistency of their directions; ||f_i|| and ||f_j|| represent the L2 norms (modulus) of the two vectors, respectively; dividing by the product of the modulus achieves cosine similarity normalization; τ is a temperature parameter that controls the smoothness of the similarity distribution, and this invention takes τ=0.1; exp(·) is the natural exponential function, which maps the cosine similarity to a positive value. When the feature vectors of two positions have the same direction (semantic similarity), the inner product is large, and the value of A(i,j) is large; when the feature vectors have different directions (semantic differences), the inner product is small or negative, and the value of A(i,j) is small. Therefore, affinity A(i,j) quantitatively describes the probability that two pixels belong to the same semantic category.

[0048] To reduce computational complexity (the full N×N affinity matrix is ​​computationally too large), this invention only calculates the affinity between pixel pairs within an 8-neighborhood in space, setting the affinity of non-adjacent pixels to 0. Row normalization is then applied to the affinity matrix to obtain the transition probability matrix: P(i,j) = A(i,j) / Σ_k A(i,k) Here, P(i,j) represents the probability of transitioning from position i to position j; Σ_k A(i,k) represents the sum of the affinities between position i and all its neighboring positions (normalization factor); after normalization, the sum of each row is 1, satisfying the requirements of the probability distribution. The transition probability matrix P describes the probability distribution of "diffusion" from one pixel to its neighboring pixels. Pixels with similar semantics have a high transition probability, while pixels with different semantics have a low transition probability, thus serving as the mathematical basis for subsequent asymptotic expansion.

[0049] The seed point set is encoded into an initial activation vector a^0∈R^N (where N is the total number of pixels in the feature map): elements at seed point positions are set to 1, and other elements are set to 0. The activation region is progressively expanded using iterative matrix multiplication to obtain the complete activation map. a^(t+1) = α × P^T × a^t + (1-α) × a^0 Where a^(t+1) represents the activation vector after the (t+1)th iteration; a^t represents the activation vector after the tth iteration; a^0 represents the initial activation vector (encoded by the seed point position); P^T represents the transpose of the transition probability matrix; α∈(0,1) is the diffusion coefficient, which controls the expansion range of each iteration, and this invention takes α=0.9; the (1-α)×a^0 term is used to maintain the connection with the initial seed point and prevent the expansion process from deviating from the original target. In each iteration, the P^T×a^t operation causes the activation value of the current activation position to be passed to the adjacent positions according to the transition probability—adjacent pixels with semantic similarity to the current activation region will obtain higher activation values, while adjacent pixels with different semantics will obtain lower activation values. After multiple iterations, the activation region gradually expands outward from the seed point, and finally stops automatically after reaching the target boundary (the semantic change at the boundary leads to a low transition probability). After T iterations (this invention takes T=15 iterations), each element value in the activation vector a^T represents the confidence that the corresponding position belongs to the bulldozing zone.

[0050] To further ensure the spatial continuity of the extended region and avoid incorrect association of semantically similar but spatially non-adjacent pixels, a spatial distance penalty is added to the affinity calculation: A_spatial(i,j) = A(i,j) × exp(-||pos_i - pos_j||² / (2r²)) Where A_spatial(i,j) represents the affinity after adding spatial constraints; A(i,j) is the original semantic affinity; pos_i and pos_j represent the spatial coordinates (two-dimensional vectors) of positions i and j, respectively; ||pos_i-pos_j|| represents the Euclidean distance between the two positions; r is the spatial bandwidth parameter, which is taken as r=3 (pixels) in this invention; exp(-d² / (2r²)) is a Gaussian decay function, with greater decay as the distance increases. The spatial distance penalty significantly reduces the affinity between pixels that are semantically similar but spatially far apart, thus ensuring the compactness and continuity of the extended region. The output of this stage is the complete region activation map a^T, which is used for boundary optimization and pseudo-label generation.

[0051] As a preferred embodiment, generating pseudo-labels corresponding to the boundary positions specifically includes: S51. Based on the preprocessed original remote sensing image, calculate the gradient magnitude of each band and fuse them to obtain the spectral gradient boundary map.

[0052] S52. Extract multi-directional and multi-scale texture features from the original remote sensing image to obtain a texture boundary map.

[0053] S53. Extract the target boundary information from the multi-scale feature map to obtain the deep feature boundary map.

[0054] S54. Fuse the spectral gradient boundary map, texture boundary map, and depth feature boundary map to generate a fused boundary map.

[0055] S55. Optimize the activation map of the complete region by fusing the boundary map to obtain the pseudo-labels corresponding to the boundary positions.

[0056] In this embodiment, refer to Figure 4 The first row of the image shows three parallel abstract boundary maps, labeled "Spectral Gradient Boundary," "Texture Boundary," and "Depth Feature Boundary," respectively. Each map uses simplified lines and regions to represent the boundary location and intensity detected by different features (represented by line thickness or color intensity). The second row shows the fusion process, where the three types of boundaries are weighted and summed (labeled with weight symbols w_s, w_t, and w_d) to obtain the fused boundary map. The third row shows the final output, where the fused boundary map is combined with the progressively expanded activation map to generate a pseudo-label mask with precise boundaries.

[0057] Based on the complete region activation map, multiple feature information is used to detect and optimize the bulldozing area boundary, generating high-quality pseudo-labels for training the segmentation network. The multi-feature boundary perception module is the third core innovation of this invention, used to refine the bulldozing area boundary. Although the complete region activation map covers the entire target area, the boundary is often relatively blurry. By fusing spectral, texture, and depth feature information from remote sensing images, the boundary is detected in multiple feature spaces and consistently fused to obtain the accurate boundary location.

[0058] The bulldozed area differs significantly from surrounding vegetation, water bodies, and other features in its spectrum, manifesting as abrupt changes in brightness and color in the imagery. Using the pre-processed original imagery, the gradient amplitude of each band is calculated and fused.

[0059] Where E_spectral represents the spectral gradient boundary map, and each pixel value represents the boundary intensity at that location; sqrt(·) represents the square root function; Σ_{c=1}^C represents the summation over all C spectral channels; I_c represents the image grayscale value of the c-th band; These represent the partial derivatives (gradients) of the image in the horizontal and vertical directions, respectively, and can be calculated using the Sobel operator. The gradient reflects the rate of change of image grayscale; drastic grayscale changes occur at boundaries, resulting in large gradient values, while uniform regions show gradual grayscale changes and small gradient values. Spectral gradient boundaries are effective for detecting boundaries of bulldozed areas and features with large spectral differences, such as vegetation and water bodies.

[0060] The surface of bulldozed areas typically exhibits a relatively uniform rough texture, which differs from the striped texture of farmland and the regular texture of buildings. A Gabor filter bank is used to extract multi-directional, multi-scale texture features, and the gradients of these texture features are calculated as texture boundaries. E_texture = || T||, where T = [G_1*I, G_2*I, ..., G_K*I] Where E_texture represents the texture boundary map; || T|| represents the gradient magnitude of the texture feature vector T; T is a K-dimensional texture feature vector; G_k (k=1,2,...,K) represents the k-th Gabor filter. This invention uses 4 directions (0°, 45°, 90°, 135°) × 3 scales = 12 filters; * represents the convolution operation; G_k*I represents the result of convolving image I with the k-th Gabor filter. Gabor filters can capture texture patterns in specific directions and scales. Through a multi-directional, multi-scale filter bank, the texture characteristics of an image can be comprehensively characterized. At the boundary where the texture changes, the gradient of the texture feature vector is larger. Texture boundary detection is more effective for detecting the boundaries of bulldozing areas and ground features with similar spectra but different textures, such as bare soil and sand.

[0061] The high-level semantic features extracted by deep neural networks contain rich target boundary information. The semantic boundary is calculated through the intermediate layer features of the multi-scale feature maps. E_deep = || F^k||, where F^k is the depth feature map of the k-th layer. Where E_deep represents the deep feature boundary map; || F^k|| represents the gradient magnitude of the k-th layer deep feature map; F^k is a multi-scale feature map, and this invention uses the conv4 layer features (rich in semantic information and with moderate resolution). Deep feature boundaries reflect the network's understanding of the target semantic category, producing strong boundary responses at locations where semantics change (i.e., the boundaries between different land cover categories). Compared to spectral and texture boundaries, deep feature boundaries are more robust to boundary detection in complex scenes and are less susceptible to interference from lighting, shadows, etc.

[0062] The above three types of boundary information are fused using learnable attention weights: E_fused = w_s × E_spectral + w_t × E_texture + w_d × E_deep Here, E_fused represents the fused boundary map; w_s, w_t, and w_d are the fusion weights for spectral, texture, and depth boundaries, respectively, satisfying w_s + w_t + w_d = 1. The weights can be learned through performance feedback on a small validation set or adaptively adjusted according to land cover type (e.g., urban areas emphasize texture boundary weights w_t, while agricultural areas emphasize spectral boundary weights w_s). Weighted fusion combines the advantages of three boundary detection methods: spectral boundaries are sensitive to color difference, texture boundaries are sensitive to texture changes, and depth boundaries are sensitive to semantic changes. The fused boundary map E_fused can more comprehensively and accurately locate the boundary positions of bulldozing areas.

[0063] The fused boundary map E_fused is used to optimize the activation map a^T of the complete region, generating the final pseudo-label Pseudo. The label generation rules are as follows: Pseudo(i,j) = 1 (positive sample), if a^T(i,j)>τ_pos and E_fused(i,j)<τ_edge Pseudo(i,j) = 0 (negative sample), if a^T(i,j) < τ_neg Pseudo(i,j) = Ignore, other cases In this diagram, Pseudo(i,j) represents the pseudo-label value at position (i,j); τ_pos is the positive sample threshold, which is the 70th percentile of the activation map values; τ_neg is the negative sample threshold, which is the 30th percentile of the activation map values; and τ_edge is the boundary threshold, which is the 50th percentile of the fused boundary map values. Only pixels with high activation values ​​(likely in the bulldozing zone) and not on the boundary (the boundary position is uncertain) are labeled as positive samples; pixels with low activation values ​​are labeled as negative samples; pixels near the boundary and those with uncertain activation values ​​in the middle range are marked as "ignore" and are not included in the loss calculation during subsequent training to avoid introducing label noise. The pseudo-label Pseudo is the final output of this stage and will be used as a supervision signal for training the semantic segmentation network.

[0064] As a preferred embodiment, the new remote sensing image is processed using a semantic segmentation network, specifically including: S61. Preprocess and extract features from the new remote sensing image to obtain a set of map tiles. Process the set of map tiles through a trained semantic segmentation network to output a probability map of the bulldozing area.

[0065] S62. After converting the bulldozing area probability map into a binary mask, the pixels in the binary mask that are greater than a preset threshold are marked as bulldozing areas, thus obtaining the bulldozing area pixel map.

[0066] S63. After performing opening and closing operations on the pixel image of the bulldozing area, remove the patches with an area smaller than a preset threshold.

[0067] S64. Convert the binary mask in the remaining patches into vector polygons, calculate the geometric properties of each patch, and obtain the bulldozing area map.

[0068] S65. Assign operational attributes to each bulldozing area patch in the bulldozing area map to obtain the bulldozing area detection results corresponding to the new remote sensing image.

[0069] As another preferred embodiment, each bulldozing patch in the bulldozing area map is assigned a business attribute, specifically including: S71. Obtain at least one first set of geographic features, wherein the first set of geographic features contains geospatial layer data with legally mandated administrative attributes.

[0070] S72. Perform spatial overlay analysis on each bulldozed area map patch and at least one set of first geographic elements, and assign at least one legal management attribute information from the set of first geographic elements to each bulldozed area map patch based on the analysis results.

[0071] S73. Based on preset rules, conduct compliance judgment on bulldozing area patches that have been assigned legal management attribute information, and generate bulldozing area detection results that include legal management attribute information and compliance judgment results.

[0072] In this embodiment, pseudo-labels are used to train a semantic segmentation network in a fully supervised manner, and inference prediction is performed on new images to output bulldozing area detection results. The semantic segmentation network adopts an encoder-decoder structure. The encoder uses a pre-trained backbone network (ResNet-50), and the parameters can be frozen or fine-tuned. The decoder adopts a lightweight design, including a feature pyramid fusion module (FPN) that upsamples and fuses multi-layer features from the encoder to restore spatial resolution; dilated spatial pyramid pooling (ASPP) to capture multi-scale contextual information and enhance the perception of targets of different sizes; and a pixel-wise classification head: a 1×1 convolutional layer that outputs a binary classification (buldozing area / background) probability map.

[0073] To address the noise inherent in weakly supervised pseudo-labels, a combined loss function is designed: L_total = L_ce + λ_1 × L_boundary + λ_2 × L_consist Where L_total is the total loss function; L_ce is the weighted cross-entropy loss; L_boundary is the boundary loss; L_consist is the consistency loss; λ_1 and λ_2 are balance coefficients, and in this invention, λ_1 = 0.5 and λ_2 = 0.3 are chosen. The weighted cross-entropy loss L_ce is defined as follows:

[0074] in, This indicates that only non-ignored pixels are summed (pixels marked "ignored" in the pseudo-labels are not included in the loss calculation); w(i,j) is the pixel weight, with a weight of 3.0 for positive samples (bulwark zone) and a weight of 1.0 for negative samples to address class imbalance; y is the pseudo-label value (0 or 1); p is the probability of the bulldozing zone predicted by the network. The boundary loss L_boundary is defined as follows: L_boundary = || pred - E_fused||² in, `pred` represents the gradient of the network's predicted probability map (larger gradients at boundary locations); `E_fused` is the fused boundary map; `||·||²` represents the square of the L2 norm. The boundary loss encourages the network's predicted boundaries to remain consistent with the multi-feature fusion boundary, improving boundary accuracy. The consistency loss `L_consist` is defined as follows: L_consist = ||pred(I) - pred(Aug(I))||² Where pred(I) represents the prediction result for the original image I; Aug(I) represents the version of image I after data augmentation (such as random flipping, rotation, color jittering); and pred(Aug(I)) represents the prediction result for the augmented image. Consistency loss requires the network to produce consistent predictions for different augmented versions of the same image, thereby improving the robustness of the model and reducing its sensitivity to false label noise.

[0075] This invention employs a two-stage training strategy. Stage 1 involves pseudo-label-guided pre-training, where the segmentation network is trained using pseudo-labels at a learning rate of 1e-4 for 50 epochs. The goal of this stage is to enable the network to learn the basic feature patterns of the bulldozing region. Stage 2 involves self-training optimization, where the network trained in Stage 1 is used to re-predict the training set. High-confidence predictions (pixels with foreground probabilities >0.9 or <0.1) are selected to update the pseudo-labels, and then training continues. This process iterates for 2-3 rounds, gradually improving the quality of the pseudo-labels and the segmentation performance.

[0076] After preprocessing and tile segmentation of the new remote sensing image to be detected, sliding window inference is performed through a trained semantic segmentation network to output a pixel-level bulldozing area probability map. The post-processing process includes: probability thresholding, converting the probability map into a binary mask (threshold set to 0.5), and marking pixels with a probability greater than the threshold as bulldozing areas; morphological operations, performing opening operations (erosion followed by dilation) to remove isolated noise points, and closing operations (dilation followed by erosion) to fill small holes inside the target; area filtering, removing patches with an area smaller than the threshold (e.g., 100 square meters) to eliminate sporadic noise interference; vector transformation, converting the raster binary mask into vector polygons, and calculating the area, perimeter, center coordinates, and circumscribed rectangle of each patch; spatial overlay analysis, performing spatial overlay with land use status maps, urban planning maps, and basic farmland protection area maps, assigning business attributes to each bulldozing area patch (e.g., location type code, whether it is within the construction land area, whether it occupies basic farmland, etc.), generating a list of suspected illegal land use, and connecting to the law enforcement inspection business system.

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

Claims

1. A method for weakly supervised detection of bulldozing areas based on high-resolution remote sensing images, characterized in that, The method includes: S11. Preprocess and extract features from the original remote sensing images to obtain a set of patches and corresponding multi-scale feature maps; S12. Perform seed point mining on the multi-scale feature map to obtain the seed point set of the bulldozing area; S13. Based on the seed point set of the bulldozing area, determine the seed points corresponding to the multi-scale feature map, and expand from the seed points to generate a complete region activation map. S14. Based on the original remote sensing image and multi-scale feature map, the activation map of the complete region is processed by the multi-feature boundary perception algorithm to obtain the pseudo label corresponding to the boundary position. S15. Based on the tile set and pseudo-labels, the new remote sensing image is processed by a semantic segmentation network to obtain the bulldozing area detection results corresponding to the new remote sensing image.

2. The method for weakly supervised detection of bulldozing areas based on high-resolution remote sensing images according to claim 1, characterized in that, Generate a set of map tiles and corresponding multi-scale feature maps, specifically including: S21. Preprocess the original remote sensing image to obtain a set of tiles, perform binary classification and labeling on the set of tiles to form a tile dataset; S22. Construct a multi-scale feature pyramid for each tile in the tile dataset and output the target image; S23. Extract features from the target image and output a multi-scale feature map.

3. The method for weakly supervised detection of bulldozing areas based on high-resolution remote sensing images according to claim 2, characterized in that, Seed point mining is performed on multi-scale feature maps, specifically including: S31. Based on the historical background category diagram, construct a background category prototype library; S32. Extract the feature vector of each spatial location in the multi-scale feature map, calculate the Euclidean distance between the feature vector and all backgrounds in the background category prototype library, and take the minimum distance as the background belonging degree of the corresponding location in the multi-scale feature map. S33. Convert the background attribution degree into a contrast significance value, and use a Gaussian function for nonlinear mapping to obtain a contrast significance map; S34. The contrast saliency maps are fused to obtain a fused saliency map. Thresholding and local maximum detection are performed on the fused saliency map to extract the seed point set of the bulldozing area.

4. The method for weakly supervised detection of bulldozing areas based on high-resolution remote sensing images according to claim 1, characterized in that, Generate a complete region activation map, specifically including: S41. Flatten the multi-scale feature map into several feature vectors, calculate the semantic similarity between any two feature vectors, and construct the semantic affinity matrix between pixels; S42. Normalize the semantic affinity matrix to obtain the transition probability matrix, encode the seed point set of the bulldozing area into the initial activation vector, and generate the activation region. S43. Based on the transition probability matrix, the activation region is expanded by iterative matrix multiplication, characterized in that a complete activation map of the region is obtained.

5. The method for weakly supervised detection of bulldozing areas based on high-resolution remote sensing images according to claim 4, characterized in that, Incorporate spatial distance penalty into the calculation of semantic similarity between any two feature vectors.

6. The method for weakly supervised detection of bulldozing areas based on high-resolution remote sensing images according to claim 1, characterized in that, Generate pseudo-labels corresponding to the boundary positions, specifically including: S51. Based on the preprocessed original remote sensing image, calculate the gradient magnitude of each band and fuse them to obtain the spectral gradient boundary map. S52. Extract multi-directional and multi-scale texture features from the original remote sensing image to obtain a texture boundary map; S53. Extract the target boundary information from the multi-scale feature map to obtain the depth feature boundary map; S54. Fuse the spectral gradient boundary map, texture boundary map, and depth feature boundary map to generate a fused boundary map; S55. Optimize the activation map of the complete region by fusing the boundary map to obtain the pseudo-labels corresponding to the boundary positions.

7. The method for weakly supervised detection of bulldozing areas based on high-resolution remote sensing images according to claim 1, characterized in that, New remote sensing images are processed using semantic segmentation networks, specifically including: S61. Preprocess and extract features from the new remote sensing image to obtain a set of map patches. Process the set of map patches through a trained semantic segmentation network to output a probability map of the bulldozing area. S62. After converting the bulldozing area probability map into a binary mask, the pixels in the binary mask that are greater than a preset threshold are marked as bulldozing areas to obtain the bulldozing area pixel map. S63. After performing opening and closing operations on the pixel image of the bulldozing area, remove the patches with an area smaller than the preset threshold. S64. Convert the binary mask in the remaining patches into vector polygons, calculate the geometric properties of each patch, and obtain the bulldozing area map; S65. Assign operational attributes to each bulldozing area patch in the bulldozing area map to obtain the bulldozing area detection results corresponding to the new remote sensing image.

8. The method for weakly supervised detection of bulldozing areas based on high-resolution remote sensing images according to claim 1, characterized in that, Assign business attributes to each bulldozing area patch in the bulldozing area map, specifically including: S71. Obtain at least one first set of geographic features, wherein the first set of geographic features contains geospatial layer data with legal management attributes; S72. Perform spatial overlay analysis on each bulldozed area map patch and at least one set of first geographic elements, and assign at least one legal management attribute information from the set of first geographic elements to each bulldozed area map patch based on the analysis results; S73. Based on preset rules, conduct compliance judgment on bulldozing area patches that have been assigned legal management attribute information, and generate bulldozing area detection results that include legal management attribute information and compliance judgment results.

Images