A method for remote sensing monitoring of construction progress based on improved shadow index

Abstract

The application discloses a kind of construction progress remote sensing monitoring methods based on improved shadow index, it is related to remote sensing image processing technical field, including obtaining research area remote sensing image sequence and carrying out time series registration, by extracting stable scene feature to carry out dynamic scene correction;Preprocessed image is converted to multiple color spaces and selects optimal channel combination;Based on local change gradient and overall gray distribution, construction area shadow extraction index is constructed;Optimal segmentation threshold is calculated to obtain shadow binary image;Extraction contour boundary point feature carries out shadow segmentation and region growth;Based on shooting parameter and building three-dimensional geometric constraint, multi-view joint optimization equation is established, and the actual height of building is obtained by iterative solution.The application can effectively handle building shielding situation, improve shadow extraction precision and building height calculation accuracy.

Description

Technical Field

[0001] This invention relates to the field of remote sensing image processing technology, and more specifically, to a remote sensing monitoring method for building construction progress based on an improved shadow index. Background Technology

[0002] Construction progress monitoring is a core component of engineering project management, and its accuracy directly impacts schedule control, resource allocation, and cost management. Traditional progress monitoring methods primarily rely on manual monitoring, but these are inefficient, error-prone, and fail to meet the monitoring needs of key investment projects. In recent years, with the rapid development of satellite remote sensing technology, non-contact engineering monitoring using satellite or aerial imagery has gradually become a research hotspot.

[0003] In construction engineering, the main building structure is a crucial part of the construction content, and the main construction process is a key aspect of monitoring construction progress. By combining remote sensing imagery with quantitative methods, the height of buildings under construction can be obtained, thereby enabling the monitoring and assessment of their construction progress. Among the quantitative methods, the shadow measurement method can extract building shadows from single-period images and combine them with shooting parameters to estimate and invert the building height, making it a mainstream method for measuring building height.

[0004] However, measuring construction progress based on shadow extraction methods still faces key challenges, namely, the accurate extraction and height measurement of building shadows. Existing shadow extraction methods (such as threshold segmentation based on HSV and C1C2C3 models, morphological processing, etc.) are mostly designed for natural scenes or urban built-up areas, making them difficult to directly apply to construction scenarios. Construction areas often have complex backgrounds and significant differences in the spectral reflectance of building materials, resulting in poor generalization ability and low edge accuracy of traditional shadow extraction methods. To address these issues, this invention proposes a remote sensing monitoring method for building construction progress based on an improved shadow index, providing a new technical means for intelligent and automated progress management of construction projects. Summary of the Invention

[0005] The purpose of this invention is to provide a remote sensing monitoring method for building construction progress based on an improved shadow index. This method solves the technical problems of inaccurate shadow extraction and large errors in building height calculation in the prior art by using multi-color space feature fusion and adaptive local and global feature weights.

[0006] This invention provides a remote sensing monitoring method for building construction progress based on an improved shading index, comprising the following steps:

[0007] The remote sensing image sequence of the study area is acquired and time-series registered. Stable scene features are extracted to construct matching point pairs. Dynamic scene correction is performed based on the matching point pairs to obtain the preprocessed image of the study area.

[0008] The preprocessed images of the study area are converted to multiple color spaces, the complementary features of each color space are calculated, the optimal channel combination is selected based on the complementary features, and the channel components are extracted.

[0009] The local variation gradient and overall gray-scale distribution of the channel components are calculated separately. The weighting coefficients are determined based on the characteristic differences between the local variation gradients and the overall gray-scale distributions. The weighting coefficients are combined with the corresponding channel components to construct the shadow extraction index of the construction area, and the shadow feature intensity map is obtained.

[0010] The optimal segmentation threshold is calculated based on the shadow feature intensity map to obtain the shadow binary image;

[0011] Extract the contour boundary points of the binary image of the shadow, calculate the direction change and curvature features of the contour boundary points, identify the shadow segmentation position, perform region growing on the segmented shadow region, and obtain independent building shadows;

[0012] The boundary contour of the building shadow is extracted based on the shooting parameters of remote sensing imagery, the three-dimensional geometric constraints of the building are constructed, and a multi-view joint optimization equation is established in combination with the shooting parameters. The actual height of the building is obtained through iterative solution, and the construction progress is calculated.

[0013] Furthermore, remote sensing image sequences of the study area are acquired and temporally registered. Stable scene features are extracted to construct matching point pairs. Dynamic scene correction is performed based on the matching point pairs to obtain preprocessed images of the study area, including:

[0014] Acquire a remote sensing image sequence of the study area, and extract the initial feature points and grayscale information from the remote sensing image sequence;

[0015] Calculate the temporal grayscale change value of the initial feature point based on the grayscale information, construct a cumulative change curve of the temporal grayscale change value in time order, analyze the change trend of the cumulative change curve, and screen out stable scene features.

[0016] Spatial clustering is performed on stable scene features, the feature point density of each cluster region is calculated, and regions with a density higher than a preset density threshold are marked as stable scene regions. Local displacement vectors are calculated based on the distribution of the stable scene regions.

[0017] Within a stable scene region, control points are selected, and the position weights and deformation weights of the control points are calculated based on the local displacement vectors. The position weights and deformation weights are then combined to construct matching point pairs.

[0018] Calculate the geometric transformation relationship between matching point pairs, establish the registration residual constraint equation between adjacent time phases, iteratively solve the registration residual constraint equation to obtain the correction parameters, and correct the remote sensing image sequence according to the correction parameters to obtain the preprocessed study area image.

[0019] Furthermore, the preprocessed images of the study area are converted to multiple color spaces, the complementary features of each color space are calculated, the optimal channel combination is selected based on the complementary features, and the channel components are extracted, including:

[0020] The preprocessed images of the study area were converted to multiple color spaces, and the channel data and grayscale information of each color space were extracted.

[0021] Based on the grayscale information, construct the grayscale distribution histogram for each channel, calculate the information entropy value of the grayscale distribution histogram, and generate complementary features containing channel data and information entropy values;

[0022] The gray-level co-occurrence matrix between channels is calculated based on complementary features. The channel correlation is calculated using the gray-level co-occurrence matrix. The channel correlation is used as a penalty term and the information entropy value is used as a gain term to construct a scoring criterion.

[0023] The channel data is combined and optimized according to the scoring criteria. The channel combination with the highest score is selected as the optimal channel combination. The feature vector of the optimal channel combination is extracted and normalized to obtain the channel components.

[0024] Furthermore, the local gradient and overall grayscale distribution of each channel component are calculated separately. Weighting coefficients are determined based on the characteristic differences between the local gradient and overall grayscale distribution. These weighting coefficients are then combined with the corresponding channel components to construct a shadow extraction index for the construction area, resulting in a shadow feature intensity map including:

[0025] Multi-directional edge detection is performed on the channel components to obtain a gradient magnitude map. The edge response value is extracted based on the gradient magnitude map to obtain the local change gradient. The gray-level frequency distribution and cumulative distribution function of the channel components are calculated to obtain the overall gray-level distribution.

[0026] A feature difference matrix is ​​constructed using the similarity between local gradient changes and overall grayscale distribution. Based on the feature difference matrix, the channel components are clustered to group similar features into shadow enhancement group and non-shadow suppression group.

[0027] A weight allocation function for the channel components in the shadow enhancement group is constructed based on the feature difference matrix. The first weight coefficient of the channel components in the shadow enhancement group is calculated through the weight allocation function. The first weight coefficient is combined with the corresponding channel component to obtain the numerator.

[0028] Based on the feature difference matrix, construct the weight constraint function of the channel components in the non-shadow suppression group, calculate the second weight coefficient of the channel components in the non-shadow suppression group, and combine the second weight coefficient with the corresponding channel components to obtain the denominator term;

[0029] The construction area shadow extraction index is constructed by dividing the numerator by the denominator, and the construction area shadow extraction index is normalized to generate a shadow feature intensity map.

[0030] Furthermore, the optimal segmentation threshold is calculated based on the shadow feature intensity map, resulting in a binary image of the shadow, including:

[0031] The shadow feature intensity map is divided into sub-blocks in an overlapping manner. The pixel mean and standard deviation of the sub-blocks are calculated. An edge response matrix is ​​constructed based on the pixel mean and standard deviation to generate a local feature description.

[0032] The feature distance matrix is ​​obtained by statistically analyzing the gray-level distribution differences of local features between adjacent sub-blocks. The main direction component of the feature distance matrix is ​​extracted to obtain the regional change features. A regional similarity function is constructed based on the regional change features. The regional similarity function is used to adaptively merge adjacent sub-blocks to generate a regional merging coefficient.

[0033] A sub-block connected graph is constructed based on the region merging coefficient. The maximum connected component of the sub-block connected graph is extracted. The maximum connected component is used as a seed region for region growth to obtain the region segmentation result.

[0034] Based on the region segmentation results, the clustering degree of pixel distribution and the gradient magnitude of the region edge are calculated for each segmented region. A baseline threshold is determined based on the clustering degree, and the baseline threshold is corrected using the gradient magnitude to obtain the optimal segmentation threshold.

[0035] The shadow feature intensity map is segmented according to the optimal segmentation threshold to obtain a binary image of the shadow.

[0036] Furthermore, the contour boundary points of the binary shadow image are extracted, the direction change and curvature features of the contour boundary points are calculated, the shadow segmentation position is identified, and region growing is performed on the segmented shadow region to obtain independent building shadows, including:

[0037] Extract the contour boundary points of the shadow binary image, obtain the coordinate sequence of the contour boundary points, calculate the displacement vector of adjacent contour boundary points in the coordinate sequence, and generate the contour boundary point direction sequence;

[0038] An overlapping sampling window is set on the direction sequence of contour boundary points, and the angle and direction accumulation value of the displacement vectors within the sampling window are calculated. The curvature matrix of the contour boundary points is constructed based on the angle and direction accumulation value.

[0039] The curvature matrix of the contour boundary points is decomposed into layers to obtain a layered curvature change map. The gradient difference and change amplitude of adjacent points in the layered curvature change map are calculated. The positions with gradient differences greater than a preset gradient difference threshold and the largest change amplitude are marked as shadow segmentation positions.

[0040] The shadow region is segmented along the shadow segmentation position, the contour shape features of the segmented shadow region are extracted, and the inter-region matching degree of the contour shape features is calculated.

[0041] Based on the inter-region matching degree, region growing is performed on the segmented shadow regions to obtain independent building shadows.

[0042] Furthermore, based on the shooting parameters of remote sensing imagery, the boundary contour of the building's shadow is extracted, and the three-dimensional geometric constraints of the building are constructed. A multi-view joint optimization equation is established in conjunction with the shooting parameters, and the actual height of the building is obtained through iterative solution. The construction progress is calculated, including:

[0043] Acquire remote sensing images and shooting parameters, decompose the remote sensing images into multiple scales according to resolution, extract the boundary response values ​​of each scale, calculate the correlation coefficient of the response values ​​of adjacent scales, use the correlation coefficient to weight and combine the boundary response values ​​to generate a boundary enhancement map, and extract the trough and peak values ​​from the boundary enhancement map as the shadow boundary contour points.

[0044] Calculate the sun projection direction based on the shooting parameters, calculate the curvature and direction of the shadow boundary contour points, construct a spatial repositioning of the contour points based on the curvature and sun projection direction, and obtain the shadow boundary contour.

[0045] Extract the directional gradient and length distribution of the shadow boundary contour, calculate the spatial mapping relationship between multi-temporal images, determine the positions of building vertices and shadow endpoints based on the spatial mapping relationship, and construct the three-dimensional geometric constraints of the building.

[0046] Based on the three-dimensional geometric constraints of the building and the shooting parameters, a multi-view joint optimization equation is established, and shadow occlusion compensation term and artifact suppression term are constructed. The actual height of the building is obtained by iterative optimization by dynamically updating the weights of the shadow occlusion compensation term and artifact suppression term.

[0047] Extract the time-series variation of the actual building height, perform piecewise fitting on the time-series variation to obtain the construction stage division results, calculate the construction rate and duration of each stage based on the division results, compare with the construction nodes in the design drawings, and calculate the building construction progress.

[0048] Furthermore, based on the building's three-dimensional geometric constraints and shooting parameters, a multi-view joint optimization equation is established, constructing shadow occlusion compensation and artifact suppression terms. The actual building height is obtained by iterative optimization through dynamically updating the weights of these terms, including:

[0049] The building vertex coordinates are obtained from the three-dimensional geometric constraints of the building. The building projection position is calculated based on the imaging angle and solar azimuth angle in the shooting parameters. A multi-view joint optimization equation is constructed based on the relationship between vertex coordinates and projection position.

[0050] The reflection intensity of the building surface is calculated using vertex coordinates. The reflection intensity is projected onto the ground to form an overlapping intensity distribution map. The occlusion connected sub-regions are determined based on the overlapping intensity distribution map, and shadow occlusion compensation terms are generated based on the occlusion connected sub-regions.

[0051] In the occluded connected sub-region, the shadow boundary is reconstructed by normal vector-guided adaptive interpolation. The reconstructed shadow boundary is projected onto the multi-view coordinate system according to the imaging angle. The positional deviation between the projected boundaries is calculated, and an artifact suppression term is generated based on the positional deviation.

[0052] Initial weights are assigned to the shadow occlusion compensation and artifact suppression terms based on the positional deviation. The weights of the compensation and suppression terms are dynamically updated based on the optimization residuals of the multi-view joint optimization equation. The actual height of the building is obtained through iterative optimization. The construction progress is then calculated based on the actual height of the building.

[0053] One technical solution provided in this embodiment of the invention is an electronic device, including: a memory, a processor, and a computer program stored in the memory and executable on the processor. When the processor executes the computer program, it implements the steps in the method described in any of the foregoing embodiments.

[0054] One technical solution provided in this embodiment of the invention is a computer-readable storage medium storing computer program instructions, which, when executed by a processor, implement the steps in the method described in any of the foregoing embodiments.

[0055] By using temporal registration and dynamic scene correction of multi-temporal images, the interference of scene changes on shadow extraction is effectively eliminated, improving the reliability of subsequent processing. The optimal channel combination is selected using complementary features of multiple color spaces to enhance the spectral discrimination of shadow areas. Adaptive weighting coefficients are constructed by combining local gradient changes and overall grayscale distribution to improve the adaptability of the shadow extraction index to complex scenes. Shadow segmentation positions are identified by the directional changes and curvature features of contour boundary points, achieving accurate segmentation of building shadows. Building height is calculated based on multi-view joint optimization equations, improving the accuracy of actual building height calculation and enabling accurate monitoring of construction progress. This method effectively solves the height calculation problem under complex conditions such as building occlusion and shadow overlap, and is suitable for construction progress monitoring in densely populated building areas. Attached Figure Description

[0056] To more clearly illustrate the technical solutions in the embodiments of the present invention 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 some embodiments of the present invention. For those skilled in the art, other drawings can be obtained based on these drawings without creative effort.

[0057] Figure 1 A flowchart of a remote sensing monitoring method for building construction progress based on an improved shadow index, provided in an embodiment of the present invention;

[0058] Figure 2 The feature calculation effect and shadow extraction result diagram based on the shadow extraction index of the construction area provided in the embodiments of the present invention;

[0059] Figure 3 A flowchart for monitoring building construction progress based on multi-temporal remote sensing provided in an embodiment of the present invention. Detailed Implementation

[0060] Exemplary embodiments will now be described in detail, examples of which are illustrated in the accompanying drawings. In the following description relating to the drawings, unless otherwise indicated, the same numerals in different drawings denote the same or similar elements.

[0061] like Figure 1 As shown, Figure 1 A remote sensing monitoring method for building construction progress based on an improved shading index is provided in this embodiment of the invention. The method includes the following steps:

[0062] The remote sensing image sequence of the study area is acquired and time-series registered. Stable scene features are extracted to construct matching point pairs. Dynamic scene correction is performed based on the matching point pairs to obtain the preprocessed image of the study area.

[0063] The preprocessed images of the study area are converted to multiple color spaces, the complementary features of each color space are calculated, the optimal channel combination is selected based on the complementary features, and the channel components are extracted.

[0064] The local variation gradient and overall grayscale distribution of the channel components are calculated separately. The weighting coefficients are determined based on the characteristic differences between the local variation gradients and the overall grayscale distribution. The weighting coefficients are combined with the corresponding channel components to construct the shadow extraction index of the construction area, and the shadow feature intensity map is obtained.

[0065] The optimal segmentation threshold is calculated based on the shadow feature intensity map to obtain the shadow binary image;

[0066] Extract the contour boundary points of the binary image of the shadow, calculate the direction change and curvature features of the contour boundary points, identify the shadow segmentation position, perform region growing on the segmented shadow region, and obtain independent building shadows;

[0067] The boundary contour of the building shadow is extracted based on the shooting parameters of remote sensing imagery, the three-dimensional geometric constraints of the building are constructed, and a multi-view joint optimization equation is established in combination with the shooting parameters. The actual height of the building is obtained through iterative solution, and the construction progress is calculated.

[0068] In one optional embodiment, a remote sensing image sequence of the study area is acquired and temporally registered, stable scene features are extracted to construct matching point pairs, and dynamic scene correction is performed based on the matching point pairs to obtain the preprocessed image of the study area, including:

[0069] Acquire a remote sensing image sequence of the study area, and extract the initial feature points and grayscale information from the remote sensing image sequence;

[0070] Calculate the temporal grayscale change value of the initial feature point based on the grayscale information, construct a cumulative change curve of the temporal grayscale change value in time order, analyze the change trend of the cumulative change curve, and screen out stable scene features.

[0071] Spatial clustering is performed on stable scene features, the feature point density of each cluster region is calculated, and regions with a density higher than a preset density threshold are marked as stable scene regions. Local displacement vectors are calculated based on the distribution of the stable scene regions.

[0072] Within a stable scene region, control points are selected, and the position weights and deformation weights of the control points are calculated based on the local displacement vectors. The position weights and deformation weights are then combined to construct matching point pairs.

[0073] Calculate the geometric transformation relationship between matching point pairs, establish the registration residual constraint equation between adjacent time phases, iteratively solve the registration residual constraint equation to obtain the correction parameters, and correct the remote sensing image sequence according to the correction parameters to obtain the preprocessed study area image.

[0074] Acquire remote sensing image sequences of the study area. These sequences typically consist of multiple images acquired at different times, which can be downloaded from remote sensing satellite databases or obtained via platforms such as drones. The image sequences should have similar or identical spatial resolutions, band combinations, and cover the same geographical area to ensure consistency in subsequent processing. The selection of the study area must consider both the research objectives and image acquisition conditions to ensure complete coverage of the target area.

[0075] After acquiring the remote sensing image sequence, initial feature points and grayscale information are extracted. The feature point extraction employs a scale-invariant feature transform algorithm, which can detect key points in the image at different scales and generate feature descriptors robust to rotation, scaling, and illumination changes. In practice, a key point detection threshold of 0.03 is set, and the number of feature points is controlled to approximately 3000 per square kilometer to ensure uniform and representative feature point distribution. For each feature point, its grayscale information is also extracted, including the grayscale value at the feature point's location and the grayscale distribution of its surrounding neighboring pixels.

[0076] Based on the extracted grayscale information, the temporal grayscale change value of the initial feature points is calculated. For each feature point, its grayscale value change throughout the entire remote sensing image sequence is analyzed. The temporal grayscale change value is obtained by calculating the difference in grayscale value of the feature point between adjacent time phases. The specific calculation formula is as follows: Where ΔGi represents the temporal grayscale change value of feature point i, Gi,t represents the grayscale value of feature point i at time t, and n represents the total number of image sequences. The calculated temporal grayscale change values ​​are used to construct a cumulative change curve in chronological order. The trend of the cumulative change curve is analyzed to identify curve segments with small slopes and stable changes; the feature points corresponding to these segments are considered stable scene features. A cumulative change threshold of 10% is set; when the cumulative change value of a feature point is below this threshold, it is marked as a stable scene feature.

[0077] Spatial clustering analysis was performed on the selected stable scene features, and a density clustering algorithm was used to group the stable scene feature points. During the clustering process, the cluster radius was set to 200 meters, and the minimum number of samples was 5. For each cluster region, its feature point density, i.e., the number of feature points per unit area, was calculated. Regions with a density higher than a preset density threshold were marked as stable scene regions. Based on the distribution of stable scene regions, a triangulation method was used to calculate local displacement vectors. These local displacement vectors reflect the relative displacement relationships between different regions, providing a basis for subsequent correction.

[0078] Control points are selected within the marked stable scene region. These control points should cover different parts of the study area and possess high feature stability. Based on the calculated local displacement vector, positional weights and deformation weights are assigned to each control point. Positional weights represent the accuracy of the control point's position during registration, while deformation weights reflect the control point's sensitivity to local deformation. Positional weights are proportional to the density of stable features around the control point, while deformation weights are proportional to the gradient of the local displacement vector change in the control point's region. The positional and deformation weights are combined to construct matching point pairs, which consist of control points in the source image and their corresponding positions in the target image.

[0079] Based on the constructed matching point pairs, the geometric transformation relationship between adjacent time phases is calculated. This geometric transformation relationship can be represented by affine transformation or polynomial transformation models. Affine transformation is suitable for smaller areas or cases with small deformations, while polynomial transformation is suitable for larger areas or cases with complex deformations. A registration residual constraint equation is established, where the residual represents the positional deviation of the matching point pair after transformation. The registration residual constraint equation is solved using an iterative optimization algorithm, such as the least squares method or robust estimation method, iteratively optimizing until the residual reaches a preset threshold or the number of iterations reaches its upper limit. After obtaining the optimal correction parameters, the remote sensing image sequence is corrected. The correction process includes geometric correction and radiometric correction to ensure that the corrected images maintain consistency in spatial location and radiometric characteristics.

[0080] This invention addresses the issue of insufficient accuracy in dynamically changing scenes using traditional registration methods by extracting stable scene features and constructing matching point pairs for dynamic scene correction. This method effectively identifies and utilizes stable scene features, improving registration accuracy and stability, reducing registration errors caused by changes in ground features, supporting efficient processing of large-scale remote sensing image sequences, and providing a reliable data foundation for subsequent change detection, target identification, and spatiotemporal analysis. It significantly enhances the reliability and accuracy of time-series analysis of remote sensing images.

[0081] In one optional embodiment, the preprocessed study area image is converted to multiple color spaces, the complementary features of each color space are calculated, the optimal channel combination is selected based on the complementary features, and the channel components are extracted, including:

[0082] The preprocessed images of the study area were converted to multiple color spaces, and the channel data and grayscale information of each color space were extracted.

[0083] Based on the grayscale information, construct the grayscale distribution histogram for each channel, calculate the information entropy value of the grayscale distribution histogram, and generate complementary features containing channel data and information entropy values;

[0084] The gray-level co-occurrence matrix between channels is calculated based on complementary features. The channel correlation is calculated using the gray-level co-occurrence matrix. The channel correlation is used as a penalty term and the information entropy value is used as a gain term to construct a scoring criterion.

[0085] The channel data is combined and optimized according to the scoring criteria. The channel combination with the highest score is selected as the optimal channel combination. The feature vector of the optimal channel combination is extracted and normalized to obtain the channel components.

[0086] Imagery of the study area is typically represented using the RGB color space. To fully extract color information from the images, it needs to be converted to other color spaces. Commonly used color spaces include HSV, HSI, Lab, and YCbCr. The HSV color space decomposes color information into three components: hue, saturation, and lightness, making it suitable for scenes with significant lighting variations. The HSI color space decomposes color information into hue, saturation, and lightness, and is not sensitive to lighting conditions. The Lab color space decomposes color information into lightness and two chromaticity components, which are independent of lightness. The YCbCr color space decomposes color information into lightness and two color difference components, making it suitable for processing video signals. In actual conversion, standard color space conversion formulas must be used to ensure conversion accuracy.

[0087] After conversion, channel data and grayscale information are extracted from each color space. Channel data refers to the component values ​​in each color space, such as the red, green, and blue channels in the RGB color space, and the hue, saturation, and lightness channels in the HSV color space. For each color space, data from all channels is extracted, and the corresponding grayscale information is calculated. Grayscale information reflects the brightness distribution of the image and can be obtained by converting the color channels into single grayscale values. During the extraction process, for high-resolution images, a region partitioning strategy can be used to divide the image into multiple sub-regions for processing, improving computational efficiency.

[0088] Based on the extracted grayscale information, a grayscale distribution histogram is constructed for each channel. The grayscale distribution histogram is an important tool for describing the grayscale distribution characteristics of an image; it statistically analyzes the frequency of occurrence of each grayscale level in the image. For an 8-bit grayscale image, the grayscale value ranges from 0 to 255, and the histogram contains 256 bars. When constructing the grayscale distribution histogram, it is necessary to count the occurrence frequency of pixels at each grayscale level and normalize it to the occurrence frequency. Based on the grayscale distribution histogram, the information entropy value of each channel is calculated. Information entropy is an indicator of the richness of image information, and its calculation formula is as follows: Here, H represents the information entropy value, L represents the number of gray levels, and ph represents the probability of a pixel with a gray value of h appearing. A higher information entropy value indicates a greater amount of information contained in the image and stronger resolution. Combining channel data and information entropy values ​​forms complementary features, which encompass both the numerical and statistical characteristics of the channels, providing a basis for subsequent channel selection.

[0089] The gray-level co-occurrence matrix (GLCM) between channels is calculated based on complementary features. The GLCM describes the spatial correlation between gray levels in an image and is an effective method for representing texture features. For any two channels, their GLCM is calculated. The element P(h, j) of the GLCM represents the frequency of occurrence of a pixel with gray value h in the first channel and a pixel with gray value j in the second channel under a specific spatial relationship. During the calculation, the spatial relationship between adjacent pixels is typically considered, including horizontal, vertical, and diagonal directions. Channel correlation is calculated using the GLCM, and this correlation can be represented by the normalized cross-correlation coefficient. In this matrix, ρuv represents the correlation coefficient between channel u and channel v, μu and μv represent the average gray values ​​of channel u and channel v, respectively, σu and σv represent the standard deviations of channel u and channel v, respectively, and P(h,j) represents the elements of the gray-level co-occurrence matrix. The closer the correlation coefficient is to 1, the stronger the correlation between the two channels.

[0090] A scoring criterion is constructed using channel correlation as a penalty and information entropy as a gain. The criterion aims to select channel combinations that are both informative and independent. For any channel combination C, its score S(C) is calculated as follows: the information entropy gain is the sum of the information entropies of all channels in the combination, and the correlation penalty is the sum of the correlation coefficients of all channel pairs in the combination. The weights of information entropy and correlation in the scoring criterion can be adjusted according to the application scenario. Generally, the information entropy weight is set to 0.7, and the correlation weight is set to 0.3 to balance the requirements of information content and independence.

[0091] The channel data is combined and optimized according to a scoring criterion, and the channel combination with the highest score is selected as the optimal channel combination. During the combination optimization process, greedy search or genetic algorithms can be used to find the highest-scoring combination among all possible channel combinations. For m color spaces with n channels, the number of possible combinations is 2n-1, which results in a large computational burden when the number of channels is large. To improve efficiency, a hierarchical search strategy can be adopted, first selecting the optimal channel within each color space, and then performing combination optimization between color spaces. After selecting the optimal channel combination, the feature vector of this combination is extracted. The feature vector contains the numerical information of each channel. The feature vector is normalized to map each channel value to a uniform range, eliminating dimensional differences and obtaining normalized channel components. Normalization can be performed using methods such as min-max normalization or Z-score standardization.

[0092] This invention achieves comprehensive capture of building shadow features by converting the study area imagery to multiple color spaces and extracting complementary features. This enhances the contrast between shadowed and unshadowed areas, improving the accuracy and reliability of construction progress monitoring. The complementary feature-based channel selection method adapts to different building environments and lighting conditions, overcoming the limitations of limited information in a single color space. It provides a high-quality data foundation for subsequent shadow extraction and construction progress analysis, significantly improving the adaptability and accuracy of remote sensing monitoring.

[0093] In one optional embodiment, the local variation gradient and overall grayscale distribution of each channel component are calculated separately. Weighting coefficients are determined based on the characteristic differences between the local variation gradients and the overall grayscale distribution. These weighting coefficients are then combined with the corresponding channel components to construct a shadow extraction index for the construction area, resulting in a shadow feature intensity map, including:

[0094] Multi-directional edge detection is performed on the channel components to obtain a gradient magnitude map. The edge response value is extracted based on the gradient magnitude map to obtain the local change gradient. The gray-level frequency distribution and cumulative distribution function of the channel components are calculated to obtain the overall gray-level distribution.

[0095] A feature difference matrix is ​​constructed using the similarity between local gradient changes and overall grayscale distribution. Based on the feature difference matrix, the channel components are clustered to group similar features into shadow enhancement group and non-shadow suppression group.

[0096] A weight allocation function for the channel components in the shadow enhancement group is constructed based on the feature difference matrix. The first weight coefficient of the channel components in the shadow enhancement group is calculated through the weight allocation function. The first weight coefficient is combined with the corresponding channel component to obtain the numerator.

[0097] Based on the feature difference matrix, construct the weight constraint function of the channel components in the non-shadow suppression group, calculate the second weight coefficient of the channel components in the non-shadow suppression group, and combine the second weight coefficient with the corresponding channel components to obtain the denominator term;

[0098] The construction area shadow extraction index is constructed by dividing the numerator by the denominator, and the construction area shadow extraction index is normalized to generate a shadow feature intensity map.

[0099] Multi-directional edge detection is performed on the channel components to obtain gradient magnitude maps. Multi-directional edge detection can employ operators such as Sobel, Prewitt, or Canny. Taking the Sobel operator as an example, convolution operations are performed on the channel components in both the horizontal and vertical directions to obtain horizontal and vertical gradient maps. The Sobel operator uses 3×3 matrices for both horizontal and vertical convolution kernels, effectively detecting edge information in the image. The gradient magnitude map is obtained by calculating the square root of the sum of the squares of the horizontal and vertical gradients. To improve detection accuracy, it can be extended to edge detection in eight directions, including angles such as 0°, 45°, 90°, and 135°. Edge response values ​​are extracted from the gradient magnitude maps. These edge response values ​​represent the edge intensity of pixels in the image; larger values ​​indicate more significant edge features. Local variation gradients are obtained by statistically analyzing the gradient values ​​of each pixel in the gradient magnitude map, reflecting the changes in the channel components within local regions.

[0100] Simultaneously calculating the gray-level frequency distribution and cumulative distribution function of each channel component yields the overall gray-level distribution. The gray-level frequency distribution statistically analyzes the frequency of each gray level in the channel components and can be plotted as a histogram. For an 8-bit quantized image, the gray level range is 0 to 255. The cumulative distribution function is the cumulative sum of the gray-level frequency distributions, reflecting the proportion of pixels less than or equal to a certain gray value among the total pixels. The formula for calculating the cumulative distribution function is: Where CDF(g) represents the cumulative probability that a gray value is less than or equal to g, and y(d) represents the probability that a pixel with a gray value of d will appear. The overall gray-level distribution combines the gray-level frequency distribution and the cumulative distribution function, comprehensively describing the gray-level characteristics of the channel components.

[0101] A feature difference matrix is ​​constructed using the similarity between local gradient changes and the overall gray-level distribution. Similarity can be calculated using methods such as correlation coefficient, mutual information, or histogram intersection. For each pair of channel components, the similarity between their local gradient changes and their overall gray-level distribution is calculated to form the feature difference matrix. The feature difference matrix is ​​a symmetric matrix, where each element represents the degree of difference between the two channel components in terms of features. Based on the feature difference matrix, the channel components are clustered using K-means clustering or hierarchical clustering algorithms. During clustering, similar features are grouped into shadow enhancement groups and non-shadow suppression groups. The shadow enhancement group contains channel components that enhance shadow features; these channel components exhibit significant gray-level differences between shadow and non-shadow regions. The non-shadow suppression group contains channel components that enhance non-shadow features; these channel components have higher response values ​​in non-shadow regions.

[0102] Based on the feature difference matrix, a weighting function is constructed for the channel components in the shadow enhancement group. The weighting function should meet the following requirements: channel components with strong responses in shadow areas should be assigned higher weights, and channel components with weak responses in shadow areas should be assigned lower weights. The weighting function can adopt an exponential decay form based on grayscale differences, and the calculation formula is as follows: In this equation, w1(f) represents the first weight coefficient of the f-th channel component in the shadow enhancement group, Ds(f) represents the grayscale difference value of the f-th channel component between the shadow and non-shadow regions, and Dmin and Dmax represent the minimum and maximum grayscale differences among all channel components, respectively. The first weight coefficient of each channel component in the shadow enhancement group is calculated using a weight allocation function, and the numerator is obtained by linearly combining the first weight coefficient with the corresponding channel component. The numerator represents the enhancement effect on shadow features; a larger value indicates a more significant shadow feature enhancement.

[0103] Based on the feature difference matrix, a weight constraint function is constructed for the channel components in the non-shadow suppression group. This function should satisfy the following requirements: assign higher weights to channel components with strong responses in non-shadow areas and lower weights to channel components with weak responses. The weight constraint function can be a normalized form based on grayscale features. The second weight coefficient of each channel component in the non-shadow suppression group is calculated, and a linear combination of the second weight coefficient and the corresponding channel component is performed to obtain the denominator. The denominator represents the enhancement effect on non-shadow features; a larger value indicates a more significant non-shadow feature. The numerator is divided by the denominator to construct the construction area shadow extraction index. This index is a ratio-based index that highlights the characteristics of shadow areas through the ratio of the numerator to the denominator. In shadow areas, a larger numerator and a smaller denominator result in a higher construction area shadow extraction index; conversely, in non-shadow areas, a smaller numerator and a larger denominator result in a lower index. To avoid calculation anomalies caused by a zero denominator, a small positive number can be added as a smoothing factor to the denominator. The shadow extraction index of the construction area is normalized, mapping the index values ​​to a range of 0 to 1, to generate a shadow feature intensity map. The normalization process can use the min-max normalization method. In the shadow feature intensity map, a pixel value closer to 1 indicates a higher probability that the point is in shadow, while a pixel value closer to 0 indicates a higher probability that the point is not in shadow. Figure 2 As shown, this embodiment demonstrates the feature calculation effect and shadow extraction result diagram based on the shadow extraction index of the construction area.

[0104] This invention constructs a shadow extraction index for construction areas by calculating the local gradient changes and overall grayscale distribution of channel components, combined with feature difference analysis and weight allocation mechanisms, thus achieving accurate shadow extraction in building construction areas. This method fully utilizes the differences in shadow characteristics across multiple color spaces, effectively distinguishing between shadowed and non-shadowed areas, and improving the accuracy and robustness of shadow extraction. Through adaptive weight adjustment, this method can adapt to shadow extraction tasks under different lighting conditions and complex backgrounds, providing reliable shadow feature data for accurate monitoring of building construction progress.

[0105] In one optional embodiment, calculating the optimal segmentation threshold based on the shadow feature intensity map to obtain the shadow binary image includes:

[0106] The shadow feature intensity map is divided into sub-blocks in an overlapping manner. The pixel mean and standard deviation of the sub-blocks are calculated. An edge response matrix is ​​constructed based on the pixel mean and standard deviation to generate a local feature description.

[0107] The feature distance matrix is ​​obtained by statistically analyzing the gray-level distribution differences of local features between adjacent sub-blocks. The main direction component of the feature distance matrix is ​​extracted to obtain the regional change features. A regional similarity function is constructed based on the regional change features. The regional similarity function is used to adaptively merge adjacent sub-blocks to generate a regional merging coefficient.

[0108] A sub-block connected graph is constructed based on the region merging coefficient. The maximum connected component of the sub-block connected graph is extracted. The maximum connected component is used as a seed region for region growth to obtain the region segmentation result.

[0109] Based on the region segmentation results, the clustering degree of pixel distribution and the gradient magnitude of the region edge are calculated for each segmented region. A baseline threshold is determined based on the clustering degree, and the baseline threshold is corrected using the gradient magnitude to obtain the optimal segmentation threshold.

[0110] The shadow feature intensity map is segmented according to the optimal segmentation threshold to obtain a binary image of the shadow.

[0111] The shadow feature intensity map is divided into sub-blocks using an overlapping sliding window strategy. Adjacent sub-blocks overlap by a certain proportion to ensure the continuity of boundary information. The sub-block size is determined based on the image resolution and shadow feature size, typically set to 32×32 pixels or 64×64 pixels. The overlap ratio is generally between 25% and 50%. A higher overlap ratio is beneficial for capturing transition features between regions but increases computational cost. For each sub-block, the pixel mean and standard deviation are calculated. The pixel mean reflects the overall brightness level of the sub-block, and the standard deviation represents the dispersion of pixel values; together, they describe the grayscale distribution characteristics of the sub-block. An edge response matrix is ​​constructed based on the pixel mean and standard deviation, representing the edge intensity distribution of pixels within the sub-block. Edge response values ​​are calculated from the grayscale difference between pixels within the sub-block and surrounding pixels; locations with large grayscale differences correspond to high edge response values, and vice versa. The edge response matrix, combined with the sub-block mean and standard deviation, generates a local feature description, which is a feature vector containing the statistical and edge features of the sub-block.

[0112] The feature distance matrix is ​​obtained by statistically analyzing the differences in gray-level distribution of local features between adjacent sub-blocks. Feature distance can be calculated using metrics such as Euclidean distance, Manhattan distance, or Mahalanobis distance. The feature distance matrix is ​​a symmetric matrix, where matrix elements represent the degree of difference in features between adjacent sub-blocks; smaller values ​​indicate greater similarity, while larger values ​​indicate less similarity. The principal direction component is extracted from the feature distance matrix to obtain the region variation features. The principal direction component can be obtained through principal component analysis, selecting the eigenvector corresponding to the largest eigenvalue of the feature distance matrix as the principal direction component. The region variation features reflect the changing trends between sub-blocks and help identify the boundaries between shaded and unshaded areas. A region similarity function is constructed based on the region variation features, and the form of the region similarity function is: Where S(q, e) represents the similarity between sub-block q and sub-block e, r(q, e) represents the feature distance between sub-block q and sub-block e, and λ is a parameter controlling the rate of change of similarity. The similarity value ranges between 0 and 1, with values ​​closer to 1 indicating greater similarity between the two sub-blocks. Adaptive merging of adjacent sub-blocks is performed using a region similarity function. When the similarity between two adjacent sub-blocks exceeds a preset threshold, they are merged into one region. The preset threshold can be dynamically adjusted according to image characteristics, typically set between 0.7 and 0.9. The merging process is iterative until no more sub-blocks can be merged. The merging result generates region merging coefficients, which describe the merging status of each sub-block with other sub-blocks.

[0113] A sub-block connected graph is constructed based on the region merging coefficient. This graph is an undirected graph where nodes represent sub-blocks and edges represent connections between adjacent sub-blocks. When two sub-blocks are merged, an edge is added to the connected graph connecting the corresponding nodes of these two sub-blocks. The maximum connected component (MPC) of the sub-block connected graph is extracted; the MPC is the subgraph with the most nodes, representing the largest homogeneous region in the image. The MPC is then used as a seed region for region growing, a segmentation method that expands progressively from the seed region. During region growing, neighboring sub-blocks similar to the seed region are added to it, with similarity determined based on a region similarity function. Region growing terminates when no more neighboring sub-blocks satisfying the similarity requirement can be added. After region growing, the resulting region segmentation divides the image into multiple homogeneous regions.

[0114] Based on the region segmentation results, the pixel distribution clustering degree and the gradient magnitude of the region edges are calculated for each segmented region. The pixel distribution clustering degree can be calculated using the variance or entropy of the pixel values ​​within the region. High clustering degree indicates uniform pixel distribution within the region, while low clustering degree indicates uneven pixel distribution. The gradient magnitude of the region edges is obtained by calculating the gradient of the pixels at the region boundaries. A large gradient magnitude indicates a clear region boundary, while a small gradient magnitude indicates a blurred region boundary. A baseline threshold is determined based on the clustering degree. The formula for the baseline threshold is: Tb = μa + α⋅σa, where Tb represents the baseline threshold, μa represents the average pixel value within the region, σa represents the standard deviation of the pixels within the region, and α is a weighting coefficient, typically between 0.5 and 2. The baseline threshold is then corrected using the gradient magnitude. The purpose of this correction is to adjust the threshold according to the clarity of the region boundaries, using a higher threshold for regions with clear boundaries and a lower threshold for regions with blurred boundaries. After correction, the optimal segmentation threshold is obtained. The optimal segmentation threshold may be different for each region, achieving adaptive segmentation.

[0115] The shadow feature intensity map is segmented based on the optimal segmentation threshold. During segmentation, pixels with intensity greater than or equal to the optimal threshold are marked as shadow regions, and pixels with intensity less than the threshold are marked as non-shadow regions. Since different thresholds may be used for each region, discontinuities in segmentation may occur at region boundaries. To address this, boundary smoothing is employed, performing local threshold interpolation at region boundaries to ensure the continuity of the segmentation results. In the resulting binary shadow image, shadow region pixels have a value of 1, and non-shadow region pixels have a value of 0, clearly representing the shadow distribution of the construction area.

[0116] This invention employs an adaptive segmentation method based on local features and region similarity to achieve accurate segmentation of shadow feature intensity maps, yielding binary shadow images with clear boundaries. This method fully considers the local characteristics and global structure of the image, overcoming the shortcomings of traditional single-threshold segmentation methods in adapting to complex scenes. It exhibits strong adaptability to uneven illumination and blurred shadow boundaries. Through region growing and boundary optimization strategies, it effectively suppresses noise and false shadow interference, improving the accuracy and robustness of shadow extraction. This provides reliable shadow region information for subsequent construction progress analysis, significantly enhancing the application effect of remote sensing monitoring in the construction field.

[0117] In one optional embodiment, the contour boundary points of the binary shadow image are extracted, the direction change and curvature features of the contour boundary points are calculated, the shadow segmentation position is identified, and region growing is performed on the segmented shadow region to obtain independent building shadows, including:

[0118] Extract the contour boundary points of the shadow binary image, obtain the coordinate sequence of the contour boundary points, calculate the displacement vector of adjacent contour boundary points in the coordinate sequence, and generate the contour boundary point direction sequence;

[0119] An overlapping sampling window is set on the direction sequence of contour boundary points, and the angle and direction accumulation value of the displacement vectors within the sampling window are calculated. The curvature matrix of the contour boundary points is constructed based on the angle and direction accumulation value.

[0120] The curvature matrix of the contour boundary points is decomposed into layers to obtain a layered curvature change map. The gradient difference and change amplitude of adjacent points in the layered curvature change map are calculated. The positions with gradient differences greater than a preset gradient difference threshold and the largest change amplitude are marked as shadow segmentation positions.

[0121] The shadow region is segmented along the shadow segmentation position, the contour shape features of the segmented shadow region are extracted, and the inter-region matching degree of the contour shape features is calculated.

[0122] Based on the inter-region matching degree, region growing is performed on the segmented shadow regions to obtain independent building shadows.

[0123] The contour boundary points of the shadow binary image are extracted using an edge tracking algorithm. Starting from the edge of the binary image, the algorithm tracks the contour point by point along the boundary between the shadow and background regions. Edge tracking begins at any edge point in the binary image and searches for the next edge point in a clockwise or counterclockwise direction until it returns to the starting point, completing the tracking of a closed contour. To ensure the continuity of the contour, an 8-neighborhood connectivity criterion is used during the search process. The coordinate sequence of the contour boundary points is obtained, arranged in the order of contour tracking to form an ordered set of points. The displacement vectors of adjacent contour boundary points in the coordinate sequence are calculated. The displacement vectors represent the direction and distance of movement from the current point to the next point. The displacement vectors are calculated by subtracting the coordinates of the previous point from the coordinates of the next point, obtaining the components in the x and y directions. A direction sequence of contour boundary points is generated based on the displacement vectors, recording the forward direction of the contour at each boundary point.

[0124] Overlapping sampling windows are set on the direction sequence of contour boundary points. The size of the sampling window is determined according to the complexity of the contour, typically set to 5% to 10% of the total contour length. The overlap rate between adjacent sampling windows is set to approximately 50% to ensure continuous contour analysis. For each sampling window, the angle and direction accumulation value of the displacement vectors within the window are calculated. The formula for calculating the angle of the displacement vectors is: Where θc represents the angle between the c-th displacement vector and the (c+1)-th displacement vector. This represents the c-th displacement vector. Let |·| represent the (c+1)th displacement vector, and |·| represent the magnitude of the displacement vector. The included angle value ranges from 0 to π; a larger included angle indicates a more abrupt turn in the contour at that point. The direction accumulation value is the sum of the directional changes of all displacement vectors within the sampling window, reflecting the degree of directional change of the contour within that region. A contour boundary point curvature matrix is ​​constructed based on the included angle and the direction accumulation value. The elements of the curvature matrix represent the curvature value of the contour at each point. A larger curvature value indicates a greater degree of curvature in the contour at that point.

[0125] The curvature matrix of the contour boundary points is decomposed hierarchically. The purpose of hierarchical decomposition is to analyze the curvature features of the contour at different scales. Hierarchical decomposition can be performed using wavelet transform or multi-scale analysis methods to decompose the curvature matrix into multiple scale levels. At lower scale levels, local detail changes in the contour can be captured; at higher scale levels, global structural features of the contour can be captured. A hierarchical curvature variation map is obtained, which shows the changes in contour curvature at different scales. The gradient difference and magnitude of change of adjacent points in the hierarchical curvature variation map are calculated. The gradient difference represents the difference in curvature values ​​between adjacent points, and the magnitude of change represents the drasticness of curvature change. The gradient difference is calculated using the difference method. For each point in the hierarchical curvature variation map, the curvature difference between it and the points before and after it is calculated. The magnitude of change can be obtained by calculating the standard deviation of the gradient difference within a local region. The location with a gradient difference greater than a preset gradient difference threshold and the largest magnitude of change is marked as the shadow segmentation location. The preset gradient difference threshold is set according to image characteristics and contour complexity, and is usually the average curvature change plus 1 to 2 times the standard deviation.

[0126] The shadow region is segmented along the segmentation line. The segmentation process starts at the segmentation point and divides the shadow region into multiple sub-regions along the segmentation line. The segmentation line can be determined using a shortest path algorithm, extending from the segmentation point into the interior of the shadow region to find the shortest path connecting different boundary points. To ensure the accuracy of the segmentation, gradient field information can be used to guide the generation of the segmentation line. The contour shape features of the segmented shadow region are extracted. These features include geometric features such as area, perimeter, compactness, and aspect ratio, as well as statistical features such as the Fourier descriptor and moment features. Contour shape feature extraction uses contour analysis methods to calculate features for each segmented shadow sub-region. The inter-region matching degree of the contour shape features is calculated, representing the degree of similarity between different shadow sub-regions. The matching degree calculation formula is: In this equation, M(Rz, Rl) represents the matching degree between regions Rz and Rl, Fk(R) represents the k-th shape feature of region R, Wk represents the weight of the k-th feature, and Q represents the total number of features. A smaller matching degree value indicates that the two regions are more similar, while a larger value indicates that the two regions are less similar.

[0127] Region growing is performed on the segmented shadow regions based on inter-region matching degree. The purpose of region growing is to merge shadow sub-regions belonging to the same building to form a complete building shadow. Region growing starts with a shadow sub-region as a seed region and gradually merges adjacent sub-regions whose matching degree with the seed region is below a threshold. The matching degree threshold is dynamically set according to the image characteristics, usually between 0.2 and 0.5. During the region growing process, the spatial relationship between shadow sub-regions must be considered, and spatially adjacent sub-regions are merged first. When there are no more sub-regions that meet the merging conditions, region growing terminates, resulting in an independent building shadow region. For the remaining unmerged shadow sub-regions, a new seed region is selected, and region growing continues until all shadow sub-regions have been processed. The final result is an independent building shadow, with each shadow region corresponding to an independent building.

[0128] This invention achieves accurate segmentation of complex shadow regions and extraction of individual building shadows by analyzing the contour boundary points of binary shadow images, combined with curvature feature recognition and region growing strategies. This method effectively handles the mutual occlusion and connection problems between building shadows, accurately identifies shadow segmentation positions, and merges related shadow regions based on contour shape features to form complete building shadows. Through multi-scale curvature analysis, this method can capture subtle changes in shadow contours, achieving high-precision shadow segmentation.

[0129] like Figure 3 As shown, in one optional embodiment, the boundary contour of the building's shadow is extracted based on the shooting parameters of remote sensing imagery, three-dimensional geometric constraints of the building are constructed, and a multi-view joint optimization equation is established in combination with the shooting parameters. The actual height of the building is obtained through iterative solution, and the construction progress is calculated, including:

[0130] Acquire remote sensing images and shooting parameters, decompose the remote sensing images into multiple scales according to resolution, extract the boundary response values ​​of each scale, calculate the correlation coefficient of the response values ​​of adjacent scales, use the correlation coefficient to weight and combine the boundary response values ​​to generate a boundary enhancement map, and extract the trough and peak values ​​from the boundary enhancement map as the shadow boundary contour points.

[0131] Calculate the sun projection direction based on the shooting parameters, calculate the curvature and direction of the shadow boundary contour points, construct a spatial repositioning of the contour points based on the curvature and sun projection direction, and obtain the shadow boundary contour.

[0132] Extract the directional gradient and length distribution of the shadow boundary contour, calculate the spatial mapping relationship between multi-temporal images, determine the positions of building vertices and shadow endpoints based on the spatial mapping relationship, and construct the three-dimensional geometric constraints of the building.

[0133] Based on the three-dimensional geometric constraints of the building and the shooting parameters, a multi-view joint optimization equation is established, and shadow occlusion compensation term and artifact suppression term are constructed. The actual height of the building is obtained by iterative optimization by dynamically updating the weights of the shadow occlusion compensation term and artifact suppression term.

[0134] Extract the time-series variation of the actual building height, perform piecewise fitting on the time-series variation to obtain the construction stage division results, calculate the construction rate and duration of each stage based on the division results, compare with the construction nodes in the design drawings, and calculate the building construction progress.

[0135] The process involves acquiring remote sensing images and acquisition parameters, including imaging time, solar altitude angle, solar azimuth angle, and satellite elevation angle. These parameters can be directly extracted from image metadata. The remote sensing images are then decomposed into multiple scales according to resolution, using Gaussian pyramids or wavelet transform methods. Gaussian pyramids generate image sequences of different resolutions by continuously applying Gaussian filtering and downsampling operations to the images. Typically, a 3- to 5-layer pyramid structure is constructed, with a resolution ratio of 2:1 between adjacent layers. Boundary response values ​​for each scale are extracted, obtained by calculating gray-level gradients or the Laplacian operator. Gray-level gradient calculation involves performing horizontal and vertical differencing operations on the image to calculate the gradient magnitude. At each scale layer, a 3×3 or 5×5 gradient operator is used to convolve the image, generating a boundary response map. The correlation coefficient between adjacent scale response values ​​is calculated; this coefficient reflects the consistency of boundary features across different scales. A higher correlation coefficient indicates better stability of the boundary feature across different scales. A boundary enhancement map is generated by weighting the boundary response values ​​using correlation coefficients, with higher correlation coefficients assigned greater weights and lower correlation coefficients assigned less weights. The trough and peak values ​​are extracted from the boundary enhancement map as the shaded boundary contour points. Peak values ​​correspond to local maxima in the boundary enhancement map, and trough values ​​correspond to local minima. The extraction process employs non-maximum suppression to ensure that the extracted boundary contour points have a certain degree of sparsity and representativeness.

[0136] The solar projection direction is calculated based on the shooting parameters. The solar projection direction is the direction in which sunlight is projected onto the ground and is directly related to the solar azimuth. The solar projection direction can be represented as the angle with true north, with clockwise being positive. The curvature and direction of the shadow boundary contour points are calculated. Curvature calculation uses a discrete curvature estimation method. For each point on the contour, several points before and after it are selected to form a local contour segment, and the curvature value is estimated using polynomial fitting or circle fitting methods. Direction calculation is based on the local tangent direction of the contour points, which can be obtained through the direction of the line connecting adjacent points or the tangent direction of the locally fitted curve. Spatial repositioning of the contour points is constructed based on curvature and solar projection direction. The purpose of spatial repositioning is to correct the position of the contour points to more accurately represent the shadow boundary. The repositioning process considers the consistency between the curvature characteristics of the contour points and the solar projection direction. Positions with large curvature usually correspond to the shadows of building corner points and require more precise positioning. After repositioning, the shadow boundary contour is obtained, which more accurately represents the boundary characteristics of building shadows.

[0137] The directional gradient and length distribution of the shadow boundary contour are extracted. The directional gradient represents the rate of change of the contour direction, and the length distribution represents the length characteristics of each contour segment. The directional gradient can be obtained by calculating the difference in the tangent direction of adjacent points on the contour. The length distribution is obtained by statistically analyzing the length information of each contour segment. The contour segment division is based on the change in directional gradient, with the location of the largest change in directional gradient serving as the boundary point of the contour segment. The spatial mapping relationship between multi-temporal images is calculated. This spatial mapping relationship describes the spatial correspondence between remote sensing images acquired at different times. The spatial mapping relationship can be calculated using feature matching and geometric transformation methods. By extracting feature points from two images, the correspondence between feature points is established, and the geometric transformation parameters are estimated. Based on the spatial mapping relationship, the positions of building vertices and shadow endpoints are determined. Building vertices are the inflection points of the building contour, and shadow endpoints are the intersection points of the shadow contour and the building contour, or the inflection points of the shadow contour. Three-dimensional geometric constraints for the building are constructed. These constraints are based on the spatial relationship between the building vertices, shadow endpoints, and the solar projection direction. The core of the geometric constraint is the relationship between the building height and the shadow length, i.e., the building height equals the shadow length multiplied by the tangent of the solar altitude angle.

[0138] A multi-view joint optimization equation is established based on the building's 3D geometric constraints and imaging parameters. Multiple views refer to images acquired at different times or angles. The joint optimization equation comprehensively considers the geometric constraints under multiple views, forming an overdetermined system of equations. Shadow occlusion compensation and artifact suppression terms are constructed. The shadow occlusion compensation term handles situations where shadows are obscured by other buildings or objects, while the artifact suppression term suppresses interference from non-building shadows. Iterative optimization is performed by dynamically updating the weights of the shadow occlusion compensation and artifact suppression terms, using gradient descent or least squares methods. The purpose of dynamically updating the weights is to gradually adjust the influence of each term during the optimization process, causing the solution to converge to the global optimum. The iteration terminates when the change in the solution is less than a preset threshold or the maximum number of iterations is reached. The actual height of the building is obtained through iterative solving; this height is the vertical distance from the ground to the top of the building.

[0139] The temporal variation of building height is extracted, referring to how building height changes over time. For continuously monitored remote sensing image sequences, the difference in building height between adjacent time points is calculated to obtain the height variation sequence. The temporal variation is then piecewise fitted to obtain the construction stage division results. Piecewise linear regression or piecewise polynomial fitting methods are used for piecewise fitting. The segmentation points are selected based on abrupt changes in the rate of change, which typically correspond to transition points between construction stages. Based on the division results, the construction rate and duration of each stage are calculated. The construction rate is the building height variation for that stage divided by the time interval, and the duration is the interval from the start to the end of that stage. The construction nodes are compared with those in the design drawings, which include key nodes such as foundation completion, main structure completion, and facade completion. The comparison method involves calculating the difference between the actual construction node time and the designed construction node time. A positive difference indicates a delay in construction progress, while a negative difference indicates an advance. Based on the comparison results, the building construction progress is calculated, expressed as the ratio of the actual completed work to the planned completed work.

[0140] This invention achieves the goal of accurately extracting the shadow boundary contour of buildings by enhancing the boundary response at multiple scales and constraining the direction of solar projection. It solves the problems of shadow occlusion and artifact interference that are difficult to handle with traditional methods by utilizing three-dimensional geometric constraints and multi-view joint optimization. Based on temporal height variation analysis, it realizes automatic division of construction stages and progress calculation. This method fully utilizes the shooting parameter information of remote sensing images and effectively improves the accuracy of building height estimation through rigorous geometric constraint modeling. The multi-view joint optimization strategy enhances the robustness and adaptability of the method. Combined with temporal analysis technology, it achieves refined monitoring of the entire construction process, providing efficient and reliable technical support for building project management.

[0141] In one optional embodiment, a multi-view joint optimization equation is established based on the building's three-dimensional geometric constraints and shooting parameters. A shadow occlusion compensation term and an artifact suppression term are constructed. The solution is obtained through iterative optimization by dynamically updating the weights of the shadow occlusion compensation term and the artifact suppression term, yielding the actual building height, including:

[0142] The building vertex coordinates are obtained from the three-dimensional geometric constraints of the building. The building projection position is calculated based on the imaging angle and solar azimuth angle in the shooting parameters. A multi-view joint optimization equation is constructed based on the relationship between vertex coordinates and projection position.

[0143] The reflection intensity of the building surface is calculated using vertex coordinates. The reflection intensity is projected onto the ground to form an overlapping intensity distribution map. The occlusion connected sub-regions are determined based on the overlapping intensity distribution map, and shadow occlusion compensation terms are generated based on the occlusion connected sub-regions.

[0144] In the occluded connected sub-region, the shadow boundary is reconstructed by normal vector-guided adaptive interpolation. The reconstructed shadow boundary is projected onto the multi-view coordinate system according to the imaging angle. The positional deviation between the projected boundaries is calculated, and an artifact suppression term is generated based on the positional deviation.

[0145] Initial weights are assigned to the shadow occlusion compensation and artifact suppression terms based on the positional deviation. The weights of the compensation and suppression terms are dynamically updated based on the optimization residuals of the multi-view joint optimization equation. The actual height of the building is obtained through iterative optimization. The construction progress is then calculated based on the actual height of the building.

[0146] The building's vertex coordinates are obtained from its 3D geometric constraints. These vertex coordinates refer to the inflection points of the building's roof outline, derived by extracting feature points from the shadow boundary outline and combining them with the solar projection relationship. Vertex coordinates are represented using a geographic coordinate system, including longitude, latitude, and elevation information. The building's projected position is calculated based on the imaging angle and solar azimuth angle from the shooting parameters. The imaging angle includes the satellite or spacecraft's pitch, roll, and yaw angles, while the solar azimuth angle represents the angle between the sun's rays and true north in the horizontal plane. The calculation of the building's projected position must consider the influence of terrain undulations and is corrected based on a digital elevation model. The projection calculation uses a ray tracing method, starting from the building's vertex and extending along the direction of the sun's rays to the ground; the intersection point is the projected position. A multi-view joint optimization equation is constructed based on the relationship between the vertex coordinates and the projected position. This equation comprehensively considers the geometric relationships of the building and its shadow from different perspectives, forming an overdetermined set of equations. The basic form of the multi-view joint optimization equation is a functional relationship between vertex height, projection distance, and solar altitude angle. A residual sum of squares minimization objective function is constructed using multi-view observation data.

[0147] The surface reflection intensity of a building is calculated using vertex coordinates. Reflection intensity is related to the building's surface material, the direction of illumination, and the viewing angle. A two-way reflection distribution function is used to calculate surface reflection intensity, considering both diffuse and specular reflection. The diffuse component is proportional to the cosine of the angle between the surface normal and the illumination direction, while the specular component is related to the geometric relationship between the surface normal, the illumination direction, and the viewing direction. The reflection intensity is projected onto the ground to form an overlapping intensity distribution map, taking into account light attenuation and occlusion effects. The overlapping intensity distribution map reflects the intensity distribution of building-reflected light received by the ground and is used to identify shadowed and potentially occluded areas. Occluded connected sub-regions are determined based on the overlapping intensity distribution map. These are areas with intensity below a threshold in the overlapping intensity distribution map that are spatially connected. The threshold is set to 30% to 50% of the normal ground reflection intensity. Region connectivity is determined using a region growing method or a connected component labeling algorithm, starting from low-intensity points and gradually merging neighboring points that meet the criteria to form complete connected sub-regions. A shadow occlusion compensation term is generated based on the occluded connected sub-regions. This compensation term corrects the height estimation error caused by shadow occlusion. The compensation term is constructed as a weighted difference between the original height estimate within the occluded area and the height estimate of the surrounding unoccluded area.

[0148] In the occluded connected sub-region, adaptive interpolation guided by normal vectors is used to reconstruct the shadow boundary, ensuring the directional continuity of the reconstructed boundary. The choice of interpolation method depends on the shape and size of the occluded region; bilinear interpolation can be used for small regions, while spline interpolation or adaptive radial basis function interpolation can be used for complex regions. During interpolation, weight allocation is related to distance and normal vector consistency, with points that are close and have consistent normal vector directions receiving higher weights. The reconstruction process uses known shadow boundary points at the edge of the occluded region as boundary conditions, obtaining the internal boundary that satisfies continuity and smoothness by solving a variational problem. The reconstructed shadow boundary is projected into a multi-view coordinate system according to the imaging angle, using perspective or affine transformation, with transformation parameters determined by the satellite or spacecraft's position, attitude, and imaging model. The positional deviation between projected boundaries is calculated, representing the spatial distance between corresponding boundary points at different viewpoints. The deviation is calculated using Euclidean or Mahalanobis distance; Euclidean distance is suitable for coordinate system normalization, while Mahalanobis distance considers the error distribution characteristics in different directions. An artifact suppression term is generated based on the positional deviation, aiming to reduce the influence of non-building shadow interference. The suppression term is constructed as a ratio function of position deviation and preset threshold; the larger the deviation, the stronger the suppression.

[0149] Initial weights are assigned to the shadow occlusion compensation and artifact suppression terms based on positional deviation. The initial weights are inversely proportional to the positional deviation; regions with smaller deviations have higher compensation weights and lower suppression weights, and vice versa. The initial weights are set between 0 and 1 to ensure the stability of the initial optimization process. The weights of the compensation and suppression terms are dynamically updated based on the optimization residuals from the multi-view joint optimization equations. The optimization residuals refer to the difference between the current estimate and the actual observation. Regions with large residuals may have severe occlusion or artifact interference, requiring corresponding weight adjustments. An adaptive strategy is used for weight updates, adjusting the weights based on the residual change trend after each iteration: increasing the compensation weight when the residual decreases and increasing the suppression weight when the residual increases or fluctuates. The actual building height is obtained through iterative optimization. The iteration process uses gradient descent or conjugate gradient methods, with a maximum number of iterations set to 50 to 100. The convergence condition is that the height change between two consecutive iterations is less than 0.5 meters or the residual change rate is less than 1%. After the iteration terminates, the solution results are robustly checked, outliers are removed, and the final building height is obtained using median filtering or weighted averaging. The construction schedule is calculated based on the actual building height, using the ratio of the current height to the design height, and corrected for in conjunction with the phased project plan.

[0150] After acquiring the actual building height, construction progress is monitored through time-series analysis. Multi-temporal remote sensing images are acquired, and the building height for each temporal phase is extracted to construct a height change time-series curve. This curve reflects the building height's trend over time and can be used to identify construction stages and calculate construction rates. Construction stage identification employs a piecewise fitting method, dividing the time-series curve into different stages based on the points of change in the height change rate. Each stage corresponds to a specific phase in the construction process, such as foundation construction, main structure construction, and decoration. The construction rate is calculated as the ratio of the height change in each stage to the time interval, expressed in meters per day. Construction progress is assessed by comparing the actual height change with the planned height change; the deviation rate represents the proportion of the difference between the actual and planned progress. A positive deviation rate indicates that the progress is ahead of schedule, while a negative rate indicates that the progress is behind schedule. Progress monitoring results are presented in graphical form, including height change curves, stage divisions, rate analysis, and progress assessments, providing decision support for project management.

[0151] This invention effectively solves the occlusion problem and artifact interference that traditional single-view methods struggle to address by combining multi-view joint optimization equations with shadow occlusion compensation and artifact suppression techniques, significantly improving the accuracy of building height estimation. Based on normal vector-guided adaptive interpolation reconstruction technology, reliable recovery of shadow boundaries in occluded areas is achieved. A dynamic weight update strategy enhances the robustness and adaptability of the optimization process. Combined with time-series analysis methods, accurate calculation of construction progress and anomaly identification are realized.

[0152] One technical solution provided in this embodiment of the invention is an electronic device, including: a memory, a processor, and a computer program stored in the memory and executable on the processor. When the processor executes the computer program, it implements the steps in the method described in any of the foregoing embodiments.

[0153] One technical solution provided in this embodiment of the invention is a computer-readable storage medium storing computer program instructions, which, when executed by a processor, implement the steps in the method described in any of the foregoing embodiments.

[0154] The specific embodiments described above are preferred embodiments of the present invention and are not intended to limit the specific scope of the present invention. The scope of the present invention includes, but is not limited to, these specific embodiments. All equivalent changes made in accordance with the shape and structure of the present invention are within the protection scope of the present invention.

Claims

1. A remote sensing monitoring method for building construction progress based on an improved shading index, characterized in that, Includes the following steps: The remote sensing image sequence of the study area is acquired and time-series registered. Stable scene features are extracted to construct matching point pairs. Dynamic scene correction is performed based on the matching point pairs to obtain the preprocessed image of the study area. The preprocessed images of the study area are converted to multiple color spaces, the complementary features of each color space are calculated, the optimal channel combination is selected based on the complementary features, and the channel components are extracted. The local variation gradient and overall gray-scale distribution of the channel components are calculated separately. The weighting coefficients are determined based on the characteristic differences between the local variation gradients and the overall gray-scale distributions. The weighting coefficients are combined with the corresponding channel components to construct the shadow extraction index of the construction area, and the shadow feature intensity map is obtained. The optimal segmentation threshold is calculated based on the shadow feature intensity map to obtain the shadow binary image; Extract the contour boundary points of the binary image of the shadow, calculate the direction change and curvature features of the contour boundary points, identify the shadow segmentation position, perform region growing on the segmented shadow region, and obtain independent building shadows; The boundary contour of the building shadow is extracted based on the shooting parameters of remote sensing imagery, the three-dimensional geometric constraints of the building are constructed, and a multi-view joint optimization equation is established in combination with the shooting parameters. The actual height of the building is obtained through iterative solution, and the construction progress is calculated.

2. The method according to claim 1, characterized in that, Remote sensing image sequences of the study area were acquired and temporally registered. Stable scene features were extracted to construct matching point pairs. Dynamic scene correction was performed based on the matching point pairs to obtain preprocessed images of the study area, including: Acquire a remote sensing image sequence of the study area, and extract the initial feature points and grayscale information from the remote sensing image sequence; Calculate the temporal grayscale change value of the initial feature point based on the grayscale information, construct a cumulative change curve of the temporal grayscale change value in time order, analyze the change trend of the cumulative change curve, and screen out stable scene features. Spatial clustering is performed on stable scene features, the feature point density of each cluster region is calculated, and regions with a density higher than a preset density threshold are marked as stable scene regions. Local displacement vectors are calculated based on the distribution of the stable scene regions. Within a stable scene region, control points are selected, and the position weights and deformation weights of the control points are calculated based on the local displacement vectors. The position weights and deformation weights are then combined to construct matching point pairs. Calculate the geometric transformation relationship between matching point pairs, establish the registration residual constraint equation between adjacent time phases, iteratively solve the registration residual constraint equation to obtain the correction parameters, and correct the remote sensing image sequence according to the correction parameters to obtain the preprocessed study area image.

3. The method according to claim 1, characterized in that, The preprocessed images of the study area were converted to multiple color spaces, and the complementary features of each color space were calculated. Based on these complementary features, the optimal channel combination was selected, and the channel components were extracted, including: The preprocessed images of the study area were converted to multiple color spaces, and the channel data and grayscale information of each color space were extracted. Based on the grayscale information, construct the grayscale distribution histogram for each channel, calculate the information entropy value of the grayscale distribution histogram, and generate complementary features containing channel data and information entropy values; The gray-level co-occurrence matrix between channels is calculated based on complementary features. The channel correlation is calculated using the gray-level co-occurrence matrix. The channel correlation is used as a penalty term and the information entropy value is used as a gain term to construct a scoring criterion. The channel data is combined and optimized according to the scoring criteria. The channel combination with the highest score is selected as the optimal channel combination. The feature vector of the optimal channel combination is extracted and normalized to obtain the channel components.

4. The method according to claim 1, characterized in that, The local gradient and overall grayscale distribution of each channel component are calculated separately. Weighting coefficients are determined based on the characteristic differences between the local gradient and overall grayscale distribution. These weighting coefficients are then combined with the corresponding channel components to construct a shadow extraction index for the construction area, resulting in a shadow feature intensity map including: Multi-directional edge detection is performed on the channel components to obtain a gradient magnitude map. The edge response value is extracted based on the gradient magnitude map to obtain the local change gradient. The gray-level frequency distribution and cumulative distribution function of the channel components are calculated to obtain the overall gray-level distribution. A feature difference matrix is ​​constructed using the similarity between local gradient changes and overall grayscale distribution. Based on the feature difference matrix, the channel components are clustered to group similar features into shadow enhancement group and non-shadow suppression group. A weight allocation function for the channel components in the shadow enhancement group is constructed based on the feature difference matrix. The first weight coefficient of the channel components in the shadow enhancement group is calculated through the weight allocation function. The first weight coefficient is combined with the corresponding channel component to obtain the numerator. Based on the feature difference matrix, construct the weight constraint function of the channel components in the non-shadow suppression group, calculate the second weight coefficient of the channel components in the non-shadow suppression group, and combine the second weight coefficient with the corresponding channel components to obtain the denominator term; The construction area shadow extraction index is constructed by dividing the numerator by the denominator, and the construction area shadow extraction index is normalized to generate a shadow feature intensity map.

5. The method according to claim 1, characterized in that, The optimal segmentation threshold is calculated based on the shadow feature intensity map, resulting in the following binary shadow image: The shadow feature intensity map is divided into sub-blocks in an overlapping manner. The pixel mean and standard deviation of the sub-blocks are calculated. An edge response matrix is ​​constructed based on the pixel mean and standard deviation to generate a local feature description. The feature distance matrix is ​​obtained by statistically analyzing the gray-level distribution differences of local features between adjacent sub-blocks. The main direction component of the feature distance matrix is ​​extracted to obtain the regional change features. A regional similarity function is constructed based on the regional change features. The regional similarity function is used to adaptively merge adjacent sub-blocks to generate a regional merging coefficient. A sub-block connected graph is constructed based on the region merging coefficient. The maximum connected component of the sub-block connected graph is extracted. The maximum connected component is used as a seed region for region growth to obtain the region segmentation result. Based on the region segmentation results, the clustering degree of pixel distribution and the gradient magnitude of the region edge are calculated for each segmented region. A baseline threshold is determined based on the clustering degree, and the baseline threshold is corrected using the gradient magnitude to obtain the optimal segmentation threshold. The shadow feature intensity map is segmented according to the optimal segmentation threshold to obtain a binary image of the shadow.

6. The method according to claim 1, characterized in that, Extract the contour boundary points of the binary shadow image, calculate the direction change and curvature features of the contour boundary points, identify the shadow segmentation position, and perform region growing on the segmented shadow region to obtain independent building shadows, including: Extract the contour boundary points of the shadow binary image, obtain the coordinate sequence of the contour boundary points, calculate the displacement vector of adjacent contour boundary points in the coordinate sequence, and generate the contour boundary point direction sequence; An overlapping sampling window is set on the direction sequence of contour boundary points, and the angle and direction accumulation value of the displacement vectors within the sampling window are calculated. The curvature matrix of the contour boundary points is constructed based on the angle and direction accumulation value. The curvature matrix of the contour boundary points is decomposed into layers to obtain a layered curvature change map. The gradient difference and change amplitude of adjacent points in the layered curvature change map are calculated. The positions with gradient differences greater than a preset gradient difference threshold and the largest change amplitude are marked as shadow segmentation positions. The shadow region is segmented along the shadow segmentation position, the contour shape features of the segmented shadow region are extracted, and the inter-region matching degree of the contour shape features is calculated. Based on the inter-region matching degree, region growing is performed on the segmented shadow regions to obtain independent building shadows.

7. The method according to claim 1, characterized in that, Based on the shooting parameters of remote sensing imagery, the boundary contour of building shadows is extracted, and the three-dimensional geometric constraints of the building are constructed. A multi-view joint optimization equation is established by combining the shooting parameters, and the actual height of the building is obtained through iterative solution. The construction progress is calculated, including: Acquire remote sensing images and shooting parameters, decompose the remote sensing images into multiple scales according to resolution, extract the boundary response values ​​of each scale, calculate the correlation coefficient of the response values ​​of adjacent scales, use the correlation coefficient to weight and combine the boundary response values ​​to generate a boundary enhancement map, and extract the trough and peak values ​​from the boundary enhancement map as the shadow boundary contour points. Calculate the sun projection direction based on the shooting parameters, calculate the curvature and direction of the shadow boundary contour points, construct a spatial repositioning of the contour points based on the curvature and sun projection direction, and obtain the shadow boundary contour. Extract the directional gradient and length distribution of the shadow boundary contour, calculate the spatial mapping relationship between multi-temporal images, determine the positions of building vertices and shadow endpoints based on the spatial mapping relationship, and construct the three-dimensional geometric constraints of the building. Based on the three-dimensional geometric constraints of the building and the shooting parameters, a multi-view joint optimization equation is established, and shadow occlusion compensation term and artifact suppression term are constructed. The actual height of the building is obtained by iterative optimization by dynamically updating the weights of the shadow occlusion compensation term and artifact suppression term. Extract the time-series variation of the actual building height, perform piecewise fitting on the time-series variation to obtain the construction stage division results, calculate the construction rate and duration of each stage based on the division results, compare with the construction nodes in the design drawings, and calculate the building construction progress.

8. The method according to claim 7, characterized in that, Based on the building's 3D geometric constraints and shooting parameters, a multi-view joint optimization equation is established. Shadow occlusion compensation and artifact suppression terms are constructed. The actual building height is obtained by iterative optimization through dynamically updating the weights of these terms. The building vertex coordinates are obtained from the three-dimensional geometric constraints of the building. The building projection position is calculated based on the imaging angle and solar azimuth angle in the shooting parameters. A multi-view joint optimization equation is constructed based on the relationship between vertex coordinates and projection position. The reflection intensity of the building surface is calculated using vertex coordinates. The reflection intensity is projected onto the ground to form an overlapping intensity distribution map. The occlusion connected sub-regions are determined based on the overlapping intensity distribution map, and shadow occlusion compensation terms are generated based on the occlusion connected sub-regions. In the occluded connected sub-region, the shadow boundary is reconstructed by normal vector-guided adaptive interpolation. The reconstructed shadow boundary is projected onto the multi-view coordinate system according to the imaging angle. The positional deviation between the projected boundaries is calculated, and an artifact suppression term is generated based on the positional deviation. Initial weights are assigned to the shadow occlusion compensation and artifact suppression terms based on the positional deviation. The weights of the compensation and suppression terms are dynamically updated based on the optimization residuals of the multi-view joint optimization equation. The actual height of the building is obtained by iterative optimization.

9. An electronic device, characterized in that, include: A memory, a processor, and a computer program stored in the memory and executable on the processor, wherein the processor, when executing the computer program, implements the steps of the method as described in any one of claims 1 to 8.

10. A computer-readable storage medium, characterized in that, The computer-readable storage medium stores computer program instructions that, when executed by a processor, implement the steps of the method as described in any one of claims 1 to 8.

Images