Method and system for visualizing road engineering progress based on unmanned aerial vehicle images

Abstract

The application provides a road engineering progress visualization method and system based on unmanned aerial vehicle images, and belongs to the technical field of remote sensing images. First, multispectral RTK unmanned aerial vehicles are used to collect images, and the collected images are corrected and divided into blocks to identify semantic candidate blocks. Through multi-scale Retinex decomposition, the illumination interference is stripped, the material spectrum and the real gradient graph are separated, the spectral graph set is constructed through weighted clustering, the joint data set excluding illumination abnormalities is generated by fusing the neighborhood gradient, the anchor point time sequence is constructed based on the data set, the stable semantic anchor points are screened through the Euclidean distance, and the affine transformation matrix of each image is iteratively solved. Based on the transformation matrix, the construction change and the illumination disturbance feature are separated, the engineering quantity is counted according to the effective pixel area, the time sequence visualization progress frame material is generated, the illumination interference can be eliminated, and the road construction time sequence progress can be accurately restored.

Description

Technical Field

[0001] This invention relates to the field of remote sensing image technology, specifically to a method and system for visualizing road construction progress based on UAV images. Background Technology

[0002] Unmanned aerial vehicle (UAV) multispectral remote sensing imagery monitoring technology has been widely applied to road construction progress control. By conducting time-series aerial surveys of construction processes such as roadbed, water-stabilized layer, and asphalt pavement, it can achieve full spatial coverage and dynamic perception of the construction status, overcoming the shortcomings of traditional manual inspections, which are characterized by low efficiency and high sampling requirements. Currently, the industry generally adopts illumination compensation models to normalize multi-temporal images to eliminate brightness differences between different aerial survey periods. This practice has been regarded as a standard step in road construction image preprocessing.

[0003] However, existing illumination compensation techniques suffer from multi-dimensional parallel defects. They uniformly attribute all temporal image differences to illumination interference and forcibly normalize pixel brightness and contrast. However, this results in lowering the high reflectance spectrum of newly laid asphalt, raising the low reflectance value of exposed subgrade, and smoothing out scattering differences caused by moisture content variations, thus compressing the inherent spectral distinguishability between different materials. Simultaneously, inconsistent gradient recalibration causes the real structural boundaries, such as construction joints and uneven compaction, to be confused with illumination residuals. The image characteristics of stable background areas are also altered due to normalization, leading to temporal drift in semantic anchor points used for registration. This results in false phenomena such as inflated construction progress, stagnation, and reverse regression, significantly reducing the accuracy of road construction progress statistics and the real reliability of dynamic visualization. Summary of the Invention

[0004] This application provides a method and system for visualizing road construction progress based on UAV images, which solves the problems mentioned in the background art.

[0005] This application provides a method and system for visualizing road construction progress based on UAV imagery.

[0006] In the first aspect, this application provides a method for visualizing road engineering progress based on UAV images, including: acquiring original images using a multispectral RTK UAV, obtaining a corrected image set and identifying potential semantic candidate blocks through image correction and road spatial stability block division;

[0007] Based on a multispectral RTK UAV and a corrected image set, global and local illumination interferences are decomposed hierarchically using a multi-scale Retinex model, and the pure material spectral reflectance component map and the true gradient component map are decoupled.

[0008] By retrieving the pure material spectral reflectance component map and the true gradient component map, and constructing an enhanced material spectral feature map set through block material spectral clustering attenuation weighted analysis, and fusing the neighborhood spectral gradient to construct a joint dataset that eliminates illumination anomalies;

[0009] Based on the joint dataset and the anchor point time series constructed by combining with the latent semantic candidate blocks, after Euclidean distance screening and cross-temporal spatial registration benchmark reconstruction, a time-stable semantic anchor point set is obtained, and the exclusive affine transformation matrix of each period image is solved by fitting affine transformation iteratively.

[0010] Based on a dedicated affine transformation matrix, a binary statistical feature separation model of construction changes and lighting disturbances is performed to obtain independent feature vectors. The effective pixel area is used as the statistical benchmark for engineering quantity to complete the visualization progress generation of complete time-series animation frame materials.

[0011] Secondly, this application provides a road engineering progress visualization system based on UAV images, including: an aerial survey preprocessing module, used to acquire original images using a multispectral RTK UAV, obtain a corrected image set and identify potential semantic candidate blocks through image correction and road spatial stability block division;

[0012] The illumination decoupling module is used to decompose global and local illumination interference through a multi-scale Retinex model based on a multispectral RTK UAV and a corrected image set, and decouple the pure material spectral reflectance component map and the real gradient component map.

[0013] The feature fusion module is used to retrieve the pure material spectral reflectance component map and the true gradient component map. Through the attenuation weighted analysis of the block material spectral clusters, it constructs an enhanced material spectral feature map set and fuses the neighborhood spectral gradient to construct a joint dataset that eliminates illumination anomalies.

[0014] The image registration module is used to construct a temporal sequence of anchor points based on the joint dataset and in combination with latent semantic candidate blocks. After Euclidean distance filtering and cross-temporal spatial registration benchmark reconstruction, a temporally stable set of semantic anchor points is obtained, and the affine transformation matrix of each period image is solved iteratively by fitting affine transformation.

[0015] The progress visualization module is used to perform binary statistical feature separation modeling of construction changes and lighting disturbances based on a dedicated affine transformation matrix, obtain independent feature vectors, and use the effective pixel area as the statistical benchmark for engineering quantity to complete the visualization progress generation of complete time-series animation frame materials.

[0016] This application has at least the following technical effects or advantages:

[0017] This invention solves the industry problem of mixed coupling between illumination disturbances and construction features in traditional road progress monitoring by employing UAV-based refined aerial survey imaging correction, dual-layer illumination decoupling, and feature purification and optimization techniques. It abandons the distortion processing method of traditional global illumination normalization, instead using a multi-scale Retinex model to layer-wise separate global and local illumination interference, preserving the original spectral and gradient features of materials. This is further combined with adaptive Gaussian weighted correction, gradient scale unification, and outlier feature removal to construct a highly robust joint dataset. Simultaneously, by using temporal feature filtering and hierarchical optimization to obtain stable semantic anchor points, and iteratively solving for a dedicated affine transformation matrix, it achieves high-precision pixel registration of multi-temporal images. This eliminates monitoring errors caused by light and shadow distortion and image misalignment from the data source and spatial reference level, significantly improving the accuracy and stability of construction feature recognition.

[0018] This solution leverages Gaussian mixture models and the expectation-maximization algorithm to achieve adaptive classification of lighting disturbances and actual construction changes, accurately removing invalid interference features. Based on pure construction features, it completes effective pixel area statistics and cumulative project quantity calculation, constructing a precise multi-process progress mapping function. By building independent rendering channels for multiple processes, it generates temporally continuous and reliable visual animation frames, effectively avoiding problems such as inflated progress, regression, and disordered layered display that are common in traditional technologies. This solution achieves quantitative statistics and dynamic visualization of the progress of multiple road construction processes, fully adapting to complex road construction scenarios and improving the objectivity and intelligence of project progress monitoring. Attached Figure Description

[0019] Figure 1 This is a flowchart illustrating the method for visualizing road construction progress based on UAV images provided in this application embodiment.

[0020] Figure 2 This is a partial schematic diagram provided in an embodiment of this application;

[0021] Figure 3 This is a line graph showing the screening of pre-stabilized anchor points provided in the embodiments of this application;

[0022] Figure 4 This is a schematic diagram of the structure of the road construction progress visualization system based on UAV images provided in this application embodiment;

[0023] The components represented by each number in the attached diagram are explained below:

[0024] Aerial survey preprocessing module 11, illumination decoupling module 12, feature fusion module 13, image registration module 14, and progress visualization module 15. Detailed Implementation

[0025] It should be noted that the terms "comprising" and "having" are intended to cover non-exclusive inclusion, for example, a process, method, system, product, or server that includes a series of steps or units is not necessarily limited to those steps or units that are explicitly listed, but may include other steps or modules that are not explicitly listed or that are inherent to such processes, methods, products, or devices.

[0026] Example 1, please refer to Figures 1 to 3 This application provides a method for visualizing the progress of road construction projects based on UAV images, including the following steps:

[0027] S100: Uses a multispectral RTK drone to acquire raw images, and obtains a corrected image set and identifies potential semantic candidate blocks through image correction and road spatial stability block division;

[0028] S200: Based on a multispectral RTK UAV and a corrected image set, global and local illumination interference are decomposed hierarchically by a multi-scale Retinex model, and the pure material spectral reflectance component map and the true gradient component map are decoupled.

[0029] S300: Retrieve the pure material spectral reflectance component map and the true gradient component map, construct an enhanced material spectral feature map set through block material spectral clustering attenuation weighted analysis, and fuse the neighborhood spectral gradient to construct a joint dataset that eliminates illumination anomalies;

[0030] S400: Based on the joint dataset and combined with the potential semantic candidate blocks to construct the anchor point time series, after Euclidean distance screening and cross-temporal spatial registration benchmark reconstruction, a time-stable semantic anchor point set is obtained, and the exclusive affine transformation matrix of each period image is solved by fitting affine transformation iteratively.

[0031] S500: Based on a dedicated affine transformation matrix, it performs binary statistical feature separation modeling of construction changes and lighting disturbances to obtain independent feature vectors. Using the effective pixel area as the statistical benchmark for engineering quantities, it completes the visualization progress generation of complete time-series animation frame materials.

[0032] In this embodiment, the complete technical chain of multispectral RTK UAV aerial survey block preprocessing, dual-layer illumination decoupling, spectral gradient feature enhancement, stable anchor point registration, and binary feature separation modeling can effectively overcome the systematic defects such as spectral convergence, gradient distortion, and registration drift caused by traditional illumination normalization processing, thereby improving the accuracy and authenticity of road engineering construction progress monitoring and visualization.

[0033] First, reliable benchmark areas are locked through image correction and stable block screening. Then, relying on the multi-scale Retinex model, global and local lighting interference is completely removed, and the pure material spectrum and real gradient features that are not affected by light and shadow are separated. For example, the inherent spectral differences of newly laid asphalt, original soil base and water-bearing construction areas can be accurately preserved, avoiding the problem of traditional algorithms smoothing out material features.

[0034] Based on this, a highly robust joint dataset is constructed by eliminating abnormal illumination samples through spectral weighting enhancement and neighborhood feature fusion. By combining temporal anchor point screening and affine transformation matrix iterative optimization, high-precision spatial alignment of multi-temporal images is achieved, completely solving the semantic drift problem of traditional registration anchor points.

[0035] Ultimately, a Gaussian mixture model is used to decouple and distinguish between construction changes and lighting disturbances. The project quantity is calculated based on the effective pixel area corresponding to pure construction features, a progress mapping function is constructed, and time-series animation frames are generated. The entire process layer by layer removes lighting interference and preserves the true features of construction entities. From data source, feature modeling, spatial registration to progress visualization, the entire chain avoids problems such as inflated progress figures, time-series oscillations, and image distortion caused by lighting distortion, significantly improving the accuracy of multi-layer road construction quantity statistics and the credibility and stability of dynamic visualization.

[0036] Specifically: S100 includes: dividing aerial survey time-series batches by road mileage, deploying encrypted flight routes with multispectral RTK drones to record the optical imaging parameter set at each shooting moment;

[0037] In practice, the design station vector file of the entire road is retrieved, the basic main flight path is generated along the longitudinal direction of the road, the heading, lateral overlap rate and standard flight altitude are uniformly set, and areas with abrupt changes in spatial morphology such as ramps, steep slopes and interchanges are marked as densified deployment targets.

[0038] Read the hardware parameters of the multispectral imaging module and RTK positioning module on the UAV, and complete the onboard power-on calibration of the sensors and positioning modules; using the starting and ending range of the densified area as the boundary, add diagonal intersecting routes between the main routes, shorten the spacing between single routes and increase the image overlap ratio, and simultaneously limit the flight altitude of the densified area to be lower than the main route altitude; verify that there are no buildings or high-voltage cables colliding or conflicting with the flight paths of the main route and the densified route, and bind each main route and densified route to the corresponding section of the stationing range;

[0039] Set a fixed shooting time interval, input the preset values ​​of solar altitude angle, exposure, and gain into the flight control shooting parameter template, and store the complete set of flight route task files after the flight route is set up. The encrypted area coordinate range and imaging parameter configuration are synchronously associated in the files for subsequent UAV automatic aerial survey execution.

[0040] Road mileage markers are linear mileage markers set along the centerline of a road, dividing road segments into sections in kilometers and meters, and are used for spatial segmentation and image indexing of the road.

[0041] Industrial drones equipped with visible light and multispectral imaging sensors and RTK differential positioning modules can simultaneously acquire image pixel radiation information and centimeter-level geographic coordinates.

[0042] In addition to the main standard flight path, the encrypted flight path is an oblique, high-density supplementary flight path added to areas with complex terrain such as ramps, slopes, and interchanges, to improve the local image overlap and reconstruction accuracy.

[0043] The optical imaging parameter set consists of equipment and environmental parameters that affect image pixel brightness and spectral values, such as solar elevation angle, exposure time, sensor gain, lens vignetting factor, and shutter parameters, which are recorded synchronously during shooting.

[0044] The original image is read and combined with the optical imaging parameter set. The radiance is solved pixel by pixel using the multispectral radiometric calibration extension algorithm. The radial distortion model is used to correct the pixel coordinate offset to obtain the corrected image set.

[0045] In practice, the sensor gain, exposure time, and vignetting correction radius parameters are read frame by frame. The drone multispectral radiometric calibration calculation is performed. The official radiometric conversion extension formula of DJI multispectral sensor is adopted. The formula is an extension of the existing commercial drone radiometric calibration model. The radiometric brightness of each visible light and multispectral channel is solved pixel by pixel to eliminate the original brightness shift caused by sensor hardware gain and exposure time.

[0046] Subsequently, a radial distortion polynomial correction model was adopted, with the image center as the origin of symmetry. Based on the distortion coefficients calibrated by the lens at the factory, the coordinate offset was corrected pixel by pixel. After correction, a distortion-free multi-channel radiometric calibration image was output. Simultaneously, the ambient illuminance and solar elevation angle values ​​corresponding to each image were extracted to construct the original illumination parameter vector of a single frame image. The radiometric calibration image, the original illumination parameter vector, and the station time-series index were bound one by one and written into the time-series repository. Only pixel-level numerical correction and parameter vector binding operations were performed. No global normalized illumination compensation for brightness and contrast was performed. All feature information of the original material spectrum and local gradient was retained, and the expansion of the corrected image set was completed.

[0047] Radiance is a physical quantity of surface-reflected radiation energy received at the sensor's entrance pupil, and it is used to unify the pixel quantization benchmark for images from different exposures and gains.

[0048] The multispectral radiometric calibration extension algorithm is a pixel-level conversion algorithm derived from the official standard linear radiometric calibration formula of DJI multispectral sensors. It introduces a lens vignetting correction term, integrates a complete set of optical element parameters such as sensor gain, exposure, and dark current, and converts the physical meaningless original sampled digital quantization values ​​of the UAV multispectral sensor into standardized radiometric brightness band by band, adapting to the multispectral multichannel synchronous imaging scenario of UAVs.

[0049] The original images are raw pixel files directly output by the sensors during the flight of the multispectral RTK drone, without any radiometric correction or geometric distortion correction.

[0050] The radial distortion model uses the principal point of the image (image center) as the origin of symmetry and constructs a polynomial through radial distortion coefficients to describe the geometric deformation law of the pixels at the edge of the lens shifting inward or outward. It is combined with bilinear interpolation to complete the inverse mapping of pixel coordinates and eliminate the inherent geometric offset error of the lens.

[0051] Each image in the corrected image set is divided into a road construction foreground area and an original landform background area using a fixed threshold coarse segmentation algorithm. The road CAD vector centerline is used as a spatial constraint to divide each image in the corrected image set into continuous blocks and align them with the centerline station number boundaries.

[0052] Each phase of the calibration image set is a single batch of complete calibration images divided according to the aerial survey time sequence; among them, the complete set of multispectral calibration images, supporting parameters, and block data obtained by a complete aerial survey of the entire road line at a fixed cycle is called a phase image.

[0053] The fixed threshold coarse segmentation algorithm is a simple segmentation algorithm that uses a fixed numerical threshold of radiance to divide pixels into two categories. It uses the average radiance of a stable background block as the threshold boundary. Pixels with radiance above the threshold are classified into one category, and those below the threshold are classified into another category. No model training is required; it only performs coarse-grained region division.

[0054] The road construction foreground area refers to the surface pixel area where the roadbed, water-stabilized soil, asphalt, slope, etc. are under construction or have been completed. It is used to limit the effective calculation range of subsequent block division, spectral feature extraction, and engineering quantity statistics, and to exclude irrelevant peripheral landforms to reduce redundant calculations.

[0055] The original landform background area refers to the original surface pixel area where no excavation, filling or paving operations have been carried out. It corresponds to the aforementioned stable background block source and is used to extract the time-invariant reference spectrum as a stable reference for illumination decoupling and cross-temporal anchor point matching.

[0056] The average radiance of multiple periods is extracted block by block, and the temporal fluctuation difference is calculated to filter blocks whose fluctuations are lower than the preset amplitude threshold for multiple consecutive periods and are not within the road construction foreground area. These blocks are then marked as potential semantic candidate blocks.

[0057] The temporal fluctuation difference refers to the value obtained by subtracting the average radiance of the same image block in two adjacent corrected images. It is the result of subtracting the average radiance of the same block in the previous period from the average radiance of the current block, and is used to quantify the radiance variation range between multiple imaging periods of the same plot.

[0058] For the difference in temporal fluctuations over multiple consecutive periods, the amplitude is taken and compared with a preset threshold. Blocks whose amplitudes are consistently below the threshold and are not within the road construction prospect area are identified as potential semantic candidate blocks. These blocks serve as the candidate base for screening temporally stable semantic anchor points.

[0059] The size of each block is uniformly set to 256×256 pixels, the block boundary is aligned with the road stationing boundary point, and the start and end stationing, slope, main road surface, and ramp area attributes of each block are recorded.

[0060] In this invention, all parameters are dimensionless by using dimensionless processing technology to remove their dimensions, and all thresholds can be obtained using the mean-standard deviation method; weight values ​​can be obtained by referring to the analytic hierarchy process.

[0061] The road CAD vector centerline is a linear vector line exported from the road design drawings. It is distributed along the center of the road and carries the station coordinates of the entire line. It is used to constrain the boundaries of image space blocks and distinguish the road construction area from the surrounding terrain.

[0062] Stable background blocks are image blocks whose original terrain has not been disturbed by excavation and spreading for a long time and whose radiation brightness has been minimal for at least two periods, serving as candidate carriers for cross-temporal registration anchor points.

[0063] In this embodiment, by adapting to the road construction scenario with a refined UAV aerial survey deployment, precise image correction, and stable block screening preprocessing scheme, the technical problems of traditional aerial photography, such as single flight path, severe image distortion, and chaotic selection of reference areas, are effectively solved, providing a high-precision and high-stability raw data base for subsequent progress monitoring.

[0064] For complex areas such as ramps and steep slopes, additional densified oblique flight paths are added, along with differentiated flight parameter settings, to make up for the shortcomings of conventional uniform flight paths in special terrain image overlap and lack of imaging details. At the same time, relying on extended multispectral radiometric calibration algorithms and radial distortion models, the deviation of equipment parameters and lens geometric distortion are accurately corrected.

[0065] For example, it can effectively correct image brightness shifts caused by differences in exposure and gain at different times, and unify the radiometric quantization standard for multi-period images.

[0066] Unlike traditional techniques that blindly perform global illumination normalization, this method only performs precise correction and parameter binding, fully preserving the original true characteristics of the pavement material's spectrum and gradient, thus avoiding the smoothing out of effective construction information. Simultaneously, it combines road stationing and CAD vector centerline constraints to complete image segmentation, accurately distinguishing the construction foreground from the original terrain background. Furthermore, based on temporal radiation fluctuations, it filters out potential semantic candidate blocks free from construction disturbances, providing a fixed and reliable stable reference benchmark for subsequent image registration and illumination decoupling. This significantly reduces computational errors caused by differences in terrain, equipment, and temporal imaging, improving the overall accuracy and stability of road construction feature extraction and progress visualization from the source.

[0067] Specifically: S200 includes: extracting the solar elevation angle corresponding to each block, dynamically configuring the Gaussian convolution kernel fusion weights, and using a multi-scale Retinex model to solve the global illumination component block by block, thereby obtaining the intermediate radiation component map and the global illumination time series matrix stripped of global illumination perturbation;

[0068] In practice, all 256×256 pixel blocks of a single-period image are retrieved. Each block is bound to the original illumination parameter vector of the current image. All blocks are traversed and the solar elevation angle values ​​stored in the vector are read one by one. A linear mapping relationship between the solar elevation angle values ​​and the weights of the multi-scale Gaussian convolution kernel is established. The higher the solar elevation angle value, the higher the weight of the large-scale Gaussian kernel. The lower the solar elevation angle value, the higher the weight of the small-scale Gaussian kernel.

[0069] For example, by setting three fixed-size Gaussian convolution kernels of 15, 45 and 75, convolution operations are performed on the radiance layer within the block to obtain three sets of independent global illumination sub-components. The three sets of sub-components are then weighted and summed according to the fusion weights obtained from the mapping to finally obtain the global illumination component map corresponding to a single block.

[0070] Specifically, a multi-scale Retinex basic model is used for component decomposition. Based on the known image imaging decomposition formula, a road image decomposition model is constructed. This model refers to the radiance obtained by calibration at any coordinate of the road image, which is equal to the result of multiplying the global large-scale ambient light component, the local small-scale terrain shadow light component, and the inherent spectral reflectance component of the surface material in sequence.

[0071] The well-known image decomposition formula refers to the brightness of the observed pixel at any coordinate of the image, which is equal to the product of the incident ambient light intensity at that location and the reflectivity of the surface object itself.

[0072] Next, pixel division is used to separate the global illumination component map from the original radiance map, resulting in an intermediate component map free from global illumination disturbances. This removes pixel brightness shifts caused by large-scale ambient illumination, generating an intermediate radiance component map for this block free from global illumination disturbances. After traversing all blocks of the current image and performing the above operations, the mean global illumination component of each block is extracted. The mean is then filled into a two-dimensional data table according to the block number and time series batch. The mean values ​​of all blocks in the current image are stored as rows in the data table. The rows of data from multiple images are superimposed sequentially to construct a global illumination time series matrix. The material reflection amplitude of the blocks is not modified or compressed throughout the process. The intermediate radiance component maps of all blocks in a single period and the row data of the mean global illumination of the current period are synchronously cached. After all time series batch operations are completed, the complete intermediate radiance component map set and the global illumination time series matrix are integrated and transmitted to the next stage.

[0073] The mean of the global illumination component represents the average ambient illumination intensity of the corresponding block in the current period. It is used to quantify the illumination level of the entire area and eliminate random fluctuations in single pixels.

[0074] Each block is pre-bound with the original illumination parameter vector of the corresponding frame. The vector contains the solar elevation angle. This value is read by traversing all image blocks and used to configure the Gaussian convolution fusion weights of this block individually, covering the illumination decomposition operation of all regions along the entire line.

[0075] The Gaussian convolution kernel fusion weight refers to the weighting coefficients assigned to the 15, 45, and 75 Gaussian convolution illumination components. These weights are dynamically adjusted based on the block solar altitude angle to adapt to the surface illumination diffusion scale under different illumination angles, accurately fit the global environmental illumination distribution, and avoid global illumination decomposition distortion caused by fixed-scale convolution.

[0076] The multi-scale Retinex model is a three-component illumination stripping model that extends the classic single-scale Retinex imaging decomposition formula. It decomposes the surface radiance into global illumination, local illumination, and material intrinsic reflection, and is used to separate large-scale environmental illumination disturbances from the standardized radiance and isolate the interference of illumination factors on the material spectral values.

[0077] The intermediate radiation component map is a pixel layer obtained by removing the global illumination component from the original radiance. It is used to eliminate the global brightness shift caused by solar altitude and overall ambient illuminance, and serves as the pre-base image for subsequent solving of local illumination residuals and extraction of pure material spectra.

[0078] The global illumination time series matrix is ​​a two-dimensional numerical matrix that stores the global illumination mean of each block as columns and aerial survey time series period as rows. It is used to record the long-term global illumination variation pattern of the entire line area. Subsequent steps can extract the illumination difference between adjacent periods to distinguish illumination fluctuations from actual construction changes.

[0079] The original radiance map is a multi-channel image layer composed of standardized radiance calculated pixel by pixel after geometric distortion correction and multispectral radiometric calibration.

[0080] Read the intermediate radiation component map corresponding to the current block to extract the average material reflection reference value of the potential semantic candidate block corresponding to the intermediate radiation component map;

[0081] By subtracting the average material reflection reference value from the median radiance of each pixel within the block, the local illumination residual of a single pixel is obtained, and a local illumination residual map is acquired. This residual represents the local illumination component caused by local shadows and minor surface undulations at the pixel, and is used to eliminate local illumination interference and solve for the pure material spectral reflectance map without illumination interference. The local illumination residuals of all pixels are combined into a two-dimensional layer, which is the local illumination residual map of the block. This residual layer is equivalent to the local illumination component.

[0082] Intermediate radiance is the pixel value obtained by dividing the original radiance pixel by the global illumination component, corresponding to the value of each point on the intermediate radiance component map; only the large-scale global ambient illumination is removed, but the local shadow illumination and surface material reflection are still retained, and the unit is the same as the original radiance and the reference material reflection value, which is used to calculate the local illumination residual.

[0083] The average material reflection benchmark value is the constant reflection mean value after the global and local illumination of the potential semantic candidate blocks are completely removed. It is obtained by averaging the pure reflection components of stable blocks over multiple periods and serves as a reference standard for surface reflection without illumination interference.

[0084] By dividing the intermediate radiance by the pixel value at the same coordinate in the local illumination residual map through pixel-by-pixel division, local illumination interference can be eliminated, and a pure material spectral reflectance component map can be obtained.

[0085] The pure material spectral reflectance component map is a pixel layer that retains only the inherent spectral reflectance characteristics of the surface material itself after removing the global and local lighting interference. It is used to completely eliminate data fluctuations caused by lighting changes and accurately distinguish different construction materials such as concrete, asphalt, and bare soil. However, it is not weighted and has cross-period lighting shift noise.

[0086] The pure material spectral reflectance component map is retrieved, and pixel gradient is solved by performing pixel gradient calculation on the pure material spectral reflectance component map. The two types of lighting layers are not used in the gradient calculation process. The gradient magnitude is not scaled or corrected during gradient calculation, and the original gradient magnitude difference is preserved to generate a true gradient component map. This component map is a two-dimensional pixel layer, which is composed of the orderly arrangement of the pixel gradient corresponding to each pixel point. It completely avoids the false edges caused by lighting changes. The gradient changes only come from the boundaries of different building materials and components, and are used for subsequent region segmentation, component boundary extraction, and DBSCAN clustering to divide construction units.

[0087] Pixel gradients are calculated using the Sobel gradient operator algorithm, representing the degree of change in the spectral reflectance value of a pixel's neighborhood. Locations of abrupt reflectance changes represent the boundaries of components, road surfaces, and construction areas, while flat areas represent the interior of homogeneous materials. Spatial contour boundaries can be located using gradients, unaffected by light and shadow.

[0088] A global illumination component map is a complete layer that combines the global illumination components of all pixels in the entire image.

[0089] Global illumination component map and local illumination residual map are used to distinguish illumination fluctuations from actual construction changes and to eliminate brightness deviations caused by illumination.

[0090] Pure material spectral reflectance component maps are used to extract the inherent spectral characteristics of different construction materials such as concrete, asphalt, and bare soil.

[0091] The true gradient component map is used to accurately identify the contours of road surfaces and components, and is unaffected by light and shadow. It is used for region division and boundary determination.

[0092] Image gradients are calculated using only pure material reflection layers that are unaffected by lighting, thus avoiding interference from changes in lighting and darkness with the boundary gradient results.

[0093] Since the information on lighting, material, and edge is coupled and mixed together, it will interfere with subsequent feature calculation, clustering, and matching. After being split into four independent layers, they can be called separately. Only the reflection map is used to calculate material features, only the lighting layer is used to analyze lighting changes, and only the gradient map is used to identify component boundaries. They do not interfere with each other and eliminate the calculation error caused by lighting.

[0094] In this embodiment, a multi-scale Retinex hierarchical decoupling algorithm that dynamically adapts to the solar elevation angle is adopted. This algorithm overcomes the shortcomings of traditional fixed convolution parameter illumination correction models, such as poor adaptability and feature aliasing. It achieves accurate separation of global and local dual-layer illumination interference and effectively solves the problem of spectral and structural feature distortion caused by light and shadow and terrain undulation in road aerial survey images.

[0095] The weights of the three levels of Gaussian convolution kernels (15, 45, and 75) are dynamically adjusted according to the solar altitude angle at different times to accurately fit the global ambient light distribution under different lighting conditions. This avoids the problem of incomplete light decomposition caused by fixed-scale convolution. For example, in the aerial survey scenario at dawn and dusk with a low solar altitude angle, the global lighting interference of the narrow shadow of the road surface can be accurately eliminated by focusing on small-scale convolution.

[0096] Meanwhile, a constant material reflection benchmark is constructed based on stable semantic candidate blocks. The residual illumination generated by local terrain and minor undulations is stripped off pixel by pixel, and the illumination component and the inherent spectral information of the material are completely separated to generate a pure material spectral reflectance component map.

[0097] Furthermore, this step calculates pixel gradients solely based on the pure material spectral layer, completely excluding the use of the illumination layer in the calculation and preserving the original gradient amplitude. This eliminates the generation of false edges caused by light and shadow, accurately restoring the true boundary gradient characteristics of different building materials such as asphalt, water-stabilized soil, and roadbed. Simultaneously, a time-series global illumination matrix is ​​constructed to fully preserve the patterns of illumination changes, enabling effective data differentiation between illumination interference and construction variations. This addresses the shortcomings of traditional techniques, such as spectral convergence and structural gradient distortion, at the root of the problem. It provides high-quality, illumination-distortion-free foundational data for subsequent feature enhancement, sample classification, and accurate progress measurement.

[0098] Specifically: S300 includes: reading the mean global illumination component of each block in the global illumination time series matrix corresponding to all aerial survey periods, performing mean calculation on the mean global illumination component to obtain the baseline mean, subtracting the baseline mean from the mean global illumination component of the current single block to obtain the deviation difference, which characterizes the degree of deviation of the current illumination from the long-term average illumination. The larger the absolute value of the difference, the stronger the illumination distortion.

[0099] Retrieve the actual gradient component map corresponding to each block. If the pixel gradient is lower than the preset gradient threshold, it means that the spectral values ​​of adjacent pixels change smoothly and there is no material switching. It is determined to be the internal region of homogeneous material. Otherwise, it means that the spectrum changes abruptly and belongs to the material boundary. It is determined to be the boundary of different building materials.

[0100] Using the deviation difference as the input variable, a Gaussian attenuation weight function is constructed using a one-dimensional Gaussian distribution. For each pixel within the homogeneous material region, the Gaussian attenuation weight function is used for weighted correction to offset the spectral value shift caused by global illumination changes in different periods, while preserving the inherent spectral differences of the material itself. This yields the corrected spectral value of a single pixel after offsetting the illumination shift and constructs an enhanced material spectral feature map. The greater the illumination deviation, the lower the assigned attenuation weight, thus weakening the spectral contribution of pixels with abnormal illumination.

[0101] Homogeneous interior refers to contiguous pixel areas covered by the same type of construction material; boundary refers to pixel bands where two or more types of materials, such as asphalt, water-stabilized, and roadbed, meet.

[0102] Multiply the original pure material spectral pixel value by the corresponding Gaussian attenuation weight to obtain the single-pixel corrected spectral value after offsetting the illumination shift.

[0103] The enhanced material spectral feature atlas is a collection of all spectral layers after full-time, full-block illumination weighted correction. The internal elements are standardized material spectral numerical layers after pixel-by-pixel illumination compensation. It relies on Gaussian attenuation weights to offset the temporal global illumination shift, weaken the interference of abnormal illumination pixels, and purify the inherent spectral differences of materials, forming a standardized spectral layer set with unified temporality and higher distinguishability of building materials. This provides a basic feature without illumination distortion for subsequent solving of pure spectral temporal differences and separating construction and illumination disturbances.

[0104] The expression for the Gaussian decay weighting function is: ,in The weights are Gaussian decay weights, with values ​​ranging from 0 to 1. This is for natural exponent calculations. The deviation difference, The standard deviation of the illumination fluctuation range refers to the range of values ​​formed by the maximum and minimum values ​​of the global illumination mean over multiple periods in a single block.

[0105] When the difference in illumination deviation is small, the weight approaches 1 infinitely, and the spectral values ​​under normal illumination are fully preserved; when the difference in deviation increases, the weight decreases smoothly and continuously, without any abrupt changes in the values, and can gently suppress the spectral distortion caused by abnormal illumination.

[0106] Based on the real gradient component map, a standardized gradient mapping layer with uniform scale and unchanged magnitude relationship is obtained by linear normalization mapping of pixel gradient intervals.

[0107] The true gradient component map in this step refers to the gradient layer calculated from the pure material spectral reflectance component map. It has not undergone spectral weighting enhancement. The gradient values ​​are generated entirely based on the original material spectrum without correction for illumination bias.

[0108] Specifically, a fixed mapping interval is defined based on the original gradient amplitude interval, and a linear normalized mapping function is set. The relative magnitude relationship of the gradient between pixels is not changed throughout the process, and the gradient values ​​are only compressed to a uniform value interval. The spectral characteristics of the enhanced material are used to distinguish different construction areas, and gradient mapping is performed separately for homogeneous material blocks to avoid cross-material gradient scale confusion. After completing the mapping conversion pixel by pixel, a standardized gradient mapping feature layer is generated.

[0109] Based on the standardized gradient mapping layer and the spectral feature atlas of the enhanced material, the mean value of the enhanced spectrum and the standardized gradient magnitude within the neighborhood are extracted and combined to form a binary joint vector, using a 3×3 neighborhood of a single pixel as the calculation unit. This vector simultaneously carries the spectral properties of the material and the contour information of the region boundary, which is used for subsequent clustering, material identification, and construction unit segmentation. A single feature cannot take into account both the material and the boundary.

[0110] The enhanced spectral mean is the average value of the spectrum after corrected illumination within a 3×3 neighborhood of a pixel. It is obtained by averaging all spectral values ​​of the 3×3 region surrounding the corresponding pixel in the spectral feature map of the enhanced material.

[0111] The standardized gradient magnitude is a boundary gradient value with a uniform scale and constant relative strength, obtained by normalizing the original true gradient through a linear fixed mapping.

[0112] The deviation difference is compared with the preset deviation threshold. If the deviation difference exceeds the preset deviation threshold, it means that the current pixel value is seriously affected by strong light or shadow, and the corresponding boundary and spectral information is distorted. The corresponding spectral gradient joint vector is then removed and does not participate in subsequent modeling. After traversing all aerial survey time series batches, an effective binary joint vector is obtained. This vector is used for global pixel clustering and region segmentation.

[0113] The global illumination time series matrix, the deviation difference of each block, and the effective binary joint vector are combined to generate a joint dataset.

[0114] In this embodiment, an integrated feature optimization mechanism combining illumination deviation quantization, Gaussian adaptive weighted correction, gradient scale normalization, and outlier feature removal effectively solves problems such as spectral numerical shifts, inconsistent gradient scales, and feature confounding failures caused by illumination temporal fluctuations in multi-temporal aerial survey images, significantly improving the reliability and distinguishability of road construction material features and boundary features.

[0115] By relying on the global illumination time-series matrix to accurately quantify the degree of illumination deviation in a single block, Gaussian attenuation weights are used to adaptively correct spectral deviations for pixels within homogeneous materials. For example, the weight of abnormal pixels in roadbed areas where strong light causes spectral distortion is automatically reduced to gently suppress illumination disturbances, fully preserve the inherent spectral differences of different construction materials, and avoid the convergence of material characteristics caused by traditional uniform correction.

[0116] Simultaneously, gradient linear normalization is achieved without altering the relative magnitude relationship of the gradients, realizing a unified scale for structural features across the entire domain and avoiding the problem of gradient scale disorder in different material blocks. By constructing a 3×3 neighborhood spectral and gradient binary joint vector, material properties and contour boundary information are simultaneously fused, and invalid feature vectors with severe illumination distortion are accurately removed, ultimately generating a high-purity and highly adaptable joint dataset.

[0117] The entire processing method balances spectral fidelity and structural uniformity, correcting numerical drift caused by temporal illumination while preserving the boundary variation characteristics of real construction materials. It completely solves the defects of traditional feature modeling, such as residual illumination noise, mixed feature dimensions, and low recognition of effective information, providing high-precision and interference-resistant standardized feature data support for subsequent pixel clustering, region segmentation, anchor point matching, and change detection.

[0118] Specifically: S400 includes: based on the joint dataset, concatenating the single-pixel corrected spectral values ​​and pixel gradients of all pixels in each latent semantic candidate block to obtain anchor binary vectors, and generating anchor time-series sequences after extracting anchor binary vectors in the same latent semantic candidate block in a time-series batch.

[0119] Based on the anchor point time series, calculate the Euclidean distance between the binary vectors of two adjacent anchor points. Traverse all time series batches and record the potential semantic candidate blocks whose Euclidean distance is lower than the preset distance threshold for three consecutive periods as pre-stable anchor points.

[0120] Based on the mean Euclidean distance of each pre-stabilized anchor point, a hierarchical screening of anchor point temporal stability is performed. Simultaneously, combined with the constraint of uniform linear spatial distribution along the road, spatial distribution optimization of anchor points is implemented, generating a set of temporally stable semantic anchor points. This provides highly reliable spatial reference control points for subsequent affine transformation matrix solving. This anchor point set is the result of a second round of fine screening of pre-stabilized anchor points. Specifically, anchor points with mean distances exceeding a threshold are removed, and anchor points are further reduced and supplemented according to road layout rules, retaining only first- and second-level qualified anchor points and eliminating inferior pre-stabilized anchor points. This step relies on road station constraints to balance the spatial distribution of anchor points, supplementing control points in weak areas such as slopes and ramps, avoiding local registration failures. Simultaneously, a standardized anchor point dataset with geographic coordinates and temporal characteristics is output as a reliable reference for solving multi-temporal affine transformation matrices, eliminating spatial misalignment caused by aerial imagery capture offsets and minor deformations, and constructing a unified registration reference across the entire domain.

[0121] Specific execution steps: Read the pre-stabilized anchor points, complete the stability classification based on the average characteristic distance of anchor points from multiple consecutive images, set three-level thresholds to remove unstable anchor points, and select first- and second-level effective stable anchor points. For example, an average Euclidean distance < 0.01 is a first-level stable anchor point, 0.01 ≤ average Euclidean distance < 0.03 is a second-level stable anchor point, and an average Euclidean distance ≥ 0.03 is removed. First-level stable anchor points are preferentially selected as the temporal alignment benchmark, and the spatial distribution of anchor points follows the linear uniform layout rule of roads.

[0122] Anchor points are evenly distributed according to station intervals, and anchor points are densified and supplemented in ramp and slope areas; the geodetic coordinates, pixel coordinates, and anchor point binary vectors of all qualified anchor points are uniformly collected, and each anchor point is bound with a unique number, station interval, and stability level; the interference of illumination distortion features is avoided throughout the process to eliminate the problem of benchmark misjudgment, and finally a complete set of temporally stable semantic anchor points is encapsulated.

[0123] The first phase of aerial survey imagery was selected as the reference imagery, and all other time-series aerial survey imagery were uniformly used as images to be registered.

[0124] Based on the temporally stable semantic anchor point set, the pixel coordinates of the corresponding stable semantic anchor points in the reference image and each image to be registered are extracted to construct coordinate pairs;

[0125] The affine transformation coordinate mapping equation is used to establish the pixel coordinate mapping relationship of the coordinate pairs. The rotation scaling coefficient and translation offset coefficient are used as the parameters to be solved. The coordinates of the stable semantic anchor point corresponding to the reference image are substituted into the affine transformation coordinate mapping equation to obtain the predicted coordinates. Combined with the pixel coordinates of the image to be registered in the corresponding coordinate pairs, the mean square of the deviation of all coordinate pairs is calculated and the square root is taken to obtain the root mean square error of the coordinates.

[0126] Using the root mean square error of coordinates as the loss function, several sets of coordinate pairs are substituted into the affine transformation coordinate mapping equation, and iterative calculations are performed until the loss function converges to a preset error threshold. The dedicated affine transformation matrix for a single image to be registered is obtained by iterative optimization. This matrix is ​​an independent transformation matrix obtained by iterative optimization of a single image to be registered. Each temporal image corresponds to the other. It is used to establish the pixel coordinate mapping relationship between the reference image and the corresponding image to be registered, correct the rotation, scaling, and translation deviations caused by shooting, and achieve cross-temporal image spatial alignment. It contains a 2-row, 2-column linear coefficient matrix and a 2-row, 1-column translation vector. The elements in the linear coefficient matrix are the rotation scaling coefficients in the X-axis and Y-axis directions, the horizontal shearing coefficients, and the vertical shearing coefficients. The elements in the translation vector are the translation offset coefficients in the X-axis and Y-axis directions.

[0127] The rotation scaling factor includes the scaling factor in the X-axis and Y-axis directions, the horizontal shearing factor, and the vertical shearing factor;

[0128] The least squares algorithm iterative optimization involves first assigning initial values ​​to the rotation, scaling, translation, and offset coefficients, then substituting them into the affine equation to calculate the predicted coordinates and solve for the root mean square error. The coefficients are continuously corrected until the error converges, and finally the optimal coefficients and transformation matrix are obtained.

[0129] Stable semantic anchors refer to qualified control points included in the temporally stable semantic anchor set after coarse screening, classification, and spatial distribution optimization. They have low spectral and gradient binary feature drift across multiple time phases, are not affected by illumination or temporary clutter, and possess long-term stable semantic features and fixed geographical locations.

[0130] The root mean square error of coordinates is used to quantify the overall deviation between the predicted pixel coordinates after affine matrix mapping and the actual anchor point pixel coordinates of the image to be registered, and to measure the registration accuracy. As a criterion for iterative optimization, the smaller the error, the better the affine matrix correction effect.

[0131] Affine transformation coordinate mapping equations are well-known basic mathematical models in the fields of computer vision and remote sensing image registration. However, this scheme abandons the conventional ordinary feature points that are easily affected by light and shadow interference, and adopts stable semantic anchor points that have been decoupled by illumination and selected and graded by time sequence as matching control points. This greatly reduces the matching deviation caused by illumination disturbance and improves the registration accuracy of multi-temporal images.

[0132] In this embodiment, semantic anchors are constructed based on the binary temporal features of spectrum and gradient. Combined with the Euclidean distance hierarchical screening and spatial balance optimization strategy, along with the iteratively optimized dedicated affine transformation matrix, the industry pain points of traditional image registration features, such as susceptibility to light and shadow interference, poor stability of control points, uneven spatial distribution, and low accuracy of temporal relative alignment, are effectively solved.

[0133] This scheme abandons the traditional matching method that relies on grayscale features. It generates anchor point time series based on the pure features after illumination decoupling. It filters pre-stable anchor points by using Euclidean distance thresholds for multiple consecutive periods and sets three levels of precision to complete fine-grained elimination. For example, drifting anchor points with a mean Euclidean distance greater than or equal to 0.03 are directly eliminated, and only highly stable first- and second-level anchor points are retained.

[0134] Simultaneously, by combining road stationing constraints with optimized anchor point distribution, and specifically increasing the density of control points in complex areas such as ramps and slopes, the problems of missing reference points and registration failure in irregular road areas encountered by traditional registration methods are solved. Furthermore, by iteratively solving the affine transformation matrix specific to each time series image using the root mean square error as the loss function, rotation, scaling, and translation deviations caused by image capture are accurately corrected, achieving pixel-level precise alignment of multiple time series images.

[0135] This method relies on stable semantic anchor points that resist light interference to establish a registration benchmark throughout the process, eliminating benchmark drift and misjudgment caused by light distortion. It constructs a unified and highly robust temporal spatial benchmark across the entire domain, providing an accurate spatial alignment basis for subsequent construction change detection and calculation of time-series differences in engineering quantities, and eliminating spatial position errors in progress monitoring from the registration level.

[0136] Because the decreased spectral discrimination of materials and gradient boundary distortion jointly compress the distribution spacing of construction changes and illumination disturbances in the feature space, the statistical distributions of the two types of signals tend to overlap, making them difficult to effectively identify using conventional classification methods. This collinearity in the feature space further propagates to the progress modeling stage, causing the progress mapping function to produce non-monotonic oscillations. This manifests as spurious changes in construction progress, such as artificial increases, stagnation, or even reverse regressions, which violate physical logic. Consequently, the accuracy of progress statistics and the reliability of dynamic visualization results are significantly reduced. Therefore, step S500 is executed.

[0137] Specifically, S500 includes: mapping the pure material spectral reflectance component map and the real gradient component map to the reference coordinate system based on a dedicated affine transformation matrix to achieve a one-to-one correspondence between pixels; solving the pure material spectral temporal difference and pixel gradient temporal difference for adjacent pixels of the same period to construct a binary difference vector; and simultaneously extracting the illumination temporal difference at the same location to obtain a ternary statistical sample; this pixel temporal ternary statistical sample is used to integrate two unrelated temporal fluctuation information, namely construction entity changes and illumination environment disturbances, into the same pixel sample to construct a three-dimensional independent statistical quantity, which serves as the complete input for fitting the bidistribution and solving the separating hyperplane of the Gaussian mixture model, so as to automatically distinguish between illumination interference and real construction changes without losing the original feature information;

[0138] The reference coordinate system is the pixel coordinate system corresponding to the reference image. All pixels of the time series layer are mapped to this coordinate system to ensure that the same ground entity has a completely consistent pixel position in each time series layer, supporting the solution of cross-time series same pixel spectrum and gradient difference.

[0139] The pure material spectral time difference is the difference between the original pure material spectral reflectance component pixel values ​​of the same pixel in two adjacent periods. It reflects the spectral changes caused by construction activities such as physical replacement and repair of road surface materials and pavement layers. The illumination has been pre-stripped and will not be affected by brightness fluctuations.

[0140] The pixel gradient temporal difference is the difference between the pixel gradients of the same pixel in two adjacent periods, representing the abrupt changes in texture caused by construction in road structure, boundaries, and smoothness; false gradients caused by illumination have been eliminated, leaving only the signal of changes in the physical shape.

[0141] The illumination temporal difference is a scalar value representing the total illumination fluctuation after superimposing the differences between two adjacent images of a single pixel. It is used to characterize the spectral numerical perturbation caused by pure illumination and serves as a priori distinguishing factor between illumination and construction-related changes. The image difference includes global illumination difference and local illumination difference.

[0142] Global illumination difference is the difference between the global illumination components of two adjacent periods for a single pixel; local illumination difference is the difference between the local illumination residuals of the same pixel in two adjacent periods for a single pixel.

[0143] A binary Gaussian mixture model was built based on ternary statistical samples to classify two types of multivariate Gaussian distributions, including samples of illumination disturbance and samples of construction changes.

[0144] Using the difference in illumination time series as a prior constraint, the expectation-maximization algorithm is used to iteratively solve for the mean vector and covariance matrix corresponding to the two types of Gaussian distributions until the bivariate Gaussian mixture model converges iteratively. After convergence, the final mean vector and covariance matrix are obtained.

[0145] The probability density function of the multivariate Gaussian distribution is used to quantify the probability of pixel time-series ternary statistical samples belonging to the Gaussian distribution of the illumination disturbance category and the construction change category respectively; after multiplying the probabilities of the two types of Gaussian distributions by their corresponding mixing weights and summing them, the overall expression of the binary Gaussian mixture model is obtained.

[0146] Based on the probability calculation results, the E-step sample posterior belonging degree estimation and the M-step mean vector and covariance matrix update iteration of the EM algorithm are completed, realizing the accurate solution of the statistical distribution parameters of the two types of changes, and providing a complete distribution basis for the subsequent derivation of the linear separating hyperplane.

[0147] The optimal linear separating hyperplane is obtained based on the final mean vector and covariance matrix. All ternary statistical samples are substituted into the optimal linear separating hyperplane to complete the binary classification, retaining the inherent statistical boundaries of the two types of changes, and generating a sample set with category labels. This sample set is a dataset formed by binding each sample with a light disturbance category label or a construction change category label after all ternary statistical samples are discriminated by the linear separating hyperplane. Because it distinguishes between the two types of disturbance samples, it provides original samples with classification labels for the subsequent generation of standardized hierarchical statistical layers.

[0148] In practice, a bicomponent Gaussian mixture model is constructed to adapt to the binary distinction scenario of construction changes and illumination disturbances. This bicomponent Gaussian mixture model contains two independent multivariate Gaussian distribution components, corresponding to the pure illumination disturbance sample category and the real construction change sample category, respectively. The illumination time series difference built into the sample is used as the prior distinction constraint. The mean vector and covariance matrix of the two Gaussian distributions are solved iteratively by relying on the expectation-maximization algorithm. After the iteration process converges and the model parameters tend to stabilize, the linear separating hyperplane that can optimally separate the two sample distributions is obtained. The hyperplane equation is composed of a linear combination of pure material spectral time series difference term, pixel gradient time series difference term and constant bias term. The pure material spectral time series difference is matched with the first weight coefficient, and the pixel gradient time series difference is matched with the second weight coefficient. The two terms are weighted and summed and then the bias coefficient is added. The overall calculation result is equal to zero.

[0149] All pixel-series ternary statistical samples are sequentially substituted into the hyperplane equation. Based on the spatial position of the sample feature landing points relative to the hyperplane, the binary categories of illumination disturbance and construction change are accurately classified. The inherent statistical distribution boundaries of the two types of change data are preserved throughout the process, avoiding the loss of feature information caused by category convergence processing. Finally, a sample set carrying accurate binary category labels is output, providing compliant and distinguishable basic sample data for subsequent binary feature separation and dataset encapsulation.

[0150] The binary Gaussian mixture model is a two-component Gaussian mixture statistical model constructed by this invention to adapt to the change characteristics of road surface time-series images. It is composed of two sets of independent multivariate Gaussian distributions corresponding to illumination disturbance and construction changes, respectively, and can accurately fit the statistical distribution law of the two types of differentiated change characteristics.

[0151] Unsupervised adaptive classification modeling is achieved by iteratively solving the distribution parameters using the EM algorithm. The linear separating hyperplane is the optimal linear segmentation interface obtained by solving the Gaussian distribution parameters of the two types of samples. With the binary difference vector as the core feature dimension, a linear discriminant relationship is constructed through weight coefficients and bias coefficients. This can accurately divide the samples of illumination disturbance and construction change in the feature space while preserving the statistical independence of the two types of changes.

[0152] The mean vector represents the central feature position of the corresponding class sample in the spectral and gradient dimensions, and the covariance matrix represents the discrete distribution range of the corresponding class sample features.

[0153] This step involves unsupervised classification of all ternary statistical samples, labeling each sample with illumination disturbance or construction change, and splitting the mixed original samples into two independent sample groups.

[0154] Based on the sample set with attached category labels, independent feature vectors corresponding to two types of multivariate Gaussian distributions are extracted. These vectors are then pruned according to category differences. The illumination disturbance feature vector retains only the illumination difference component; the construction change feature vector retains both the material spectrum and gradient difference components, respectively. Irrelevant feature dimensions are removed to isolate false illumination fluctuations from road surface modification signals, providing a clean classification data source for the hierarchical construction of construction-specific feature maps and illumination interference statistical layers.

[0155] Because the independent feature vectors corresponding to the construction change samples completely isolate false pixel signals caused by light and darkness, they only carry the real construction change information of the road surface layer. By comparing the feature intervals of the process, pixels can be screened, which can avoid the over-counting or under-counting of effective area caused by light distortion, and ensure that the statistical benchmark of the project quantity is true and reliable. Other features containing light components do not participate in pixel judgment.

[0156] Extract the independent feature vectors corresponding to the construction change samples, and use the effective pixel area within the pixel region that conforms to the feature interval of the corresponding process as the statistical benchmark for the engineering quantity. This is to convert the construction coverage pixel area on the image into the actual road construction paving area. This benchmark is a unified measurement standard for the cumulative calculation of all engineering quantities, ensuring that the calculation caliber of engineering quantities is consistent for different time sequences and different processes.

[0157] A pixel region is a single-pixel geographic unit that is divided from the road surface after aerial image registration; the pixel region corresponding to a process is the set of all pixels within the road station range in a single time-series image that belong to the construction scope of that process in geographic space.

[0158] The process refers to the three-layer construction process of the road, which is in the following order: roadbed construction process, water-stabilized base course construction process, and asphalt surface course construction process.

[0159] The process feature range is a predefined numerical range, determined by the distribution range of material spectral time difference and native fidelity gradient time difference corresponding to the standard construction of each process; it is pre-calibrated based on the statistical distribution of pure construction feature vectors of the previous standard construction samples, and is used to distinguish which construction process a pixel belongs to.

[0160] The effective pixel area is calculated by iterating through each pixel within the process's pixel region and comparing its pure construction feature vector to see if it falls within the preset feature range for that process. Pixels that successfully match are considered effective pixels. The actual geographical area of ​​each single pixel is then overlaid to obtain the effective pixel area for a single time-series batch. This effective pixel area accurately corresponds to the actual construction area, avoiding false pixels caused by lighting interference and objectively reflecting the actual construction volume.

[0161] The pixel region corresponding to the process is traversed in each time sequence batch. The indicator function is used to determine whether a single pixel matches the feature interval of the corresponding process. If the match is successful, it is included in the effective construction area. If the match fails, it is not counted. The effective pixel areas of all collection periods before the current time sequence batch are accumulated in sequence and the area integration calculation is completed to obtain the cumulative engineering volume of the corresponding process up to the current time sequence batch.

[0162] The cumulative project volume represents the total construction and paving volume of the corresponding process from the start of construction to the current data collection period. It serves as the core input base for the progress mapping function and is used to convert standardized comprehensive construction progress values.

[0163] Based on the cumulative workload, the current completion percentage of a single process is obtained. After weighted summation of multiple processes, a progress mapping function is obtained. This function comprehensively reflects the overall construction progress of the section after the superposition of multiple processes such as roadbed, water-stabilized base course, and asphalt surface course. It avoids the problems of false progress inflation and regression distortion caused by light fluctuations. As the driving basis for the rendering of the visualization layers of subsequent layered processes, it controls the activation of color blocks and progress display effects of each construction layer, and supports the generation of time-series dynamic construction animation frames.

[0164] The current completion percentage of a single process is obtained by dividing the cumulative amount of work done in the single process by the total design amount of work done in that process. The total design amount of work done in the process is obtained based on the road construction design drawings, the paving width of each section at the station number, and the design paving thickness.

[0165] Based on the progress mapping function, using the pure material spectral reflectance component map and the real gradient component map as the visualization base, an independent rendering channel and layered progress color blocks are constructed for each process. The corresponding pixel color blocks are activated based on the cumulative engineering quantity judgment condition. Only the output value of the construction progress mapping function is used to control the layer rendering. This step is used to build an independent visualization channel for multiple processes and generate layered rendering layer materials with no lighting false color blocks.

[0166] Independent rendering channels, layered progress color blocks, geographic anchor points, and station numbers are loaded in batches to generate visual animation frames, resulting in complete time-series animation frame materials to visualize the dynamic changes in the construction progress of multiple road construction processes throughout the entire time sequence.

[0167] In this embodiment, a ternary statistical sample of spectrum, gradient, and illumination is constructed based on the precisely registered temporal images. By combining the Gaussian mixture model and the EM iterative algorithm, unsupervised and accurate separation of illumination disturbances and construction changes is achieved. This solves the problems of traditional technologies being unable to distinguish between light and shadow fluctuations and actual construction changes, as well as the distortion of engineering quantity statistics and the oscillation distortion of progress visualization.

[0168] By solving for the optimal linear separating hyperplane, binary classification of samples is completed, and two independent feature vectors are extracted in a differentiated manner. The dimension of light interference is completely removed, and only the spectrum and gradient changes are retained to represent the real construction state. For example, it can accurately distinguish the shadow fluctuations of asphalt surface layer from the pixel feature changes brought about by actual paving construction, and eliminate the problem of over-calculation and under-calculation of engineering quantities caused by misjudgment of light.

[0169] Meanwhile, the feature range of each process is matched with the pure construction feature vector, the effective pixel area is selected and the cumulative project volume is calculated by integration. A high-precision progress mapping function is constructed by weighting multiple processes, which avoids the defects of traditional progress models that cause progress regression and false increases due to lighting noise.

[0170] Ultimately, relying on a dedicated rendering channel and layered progress color blocks to generate time-series animation frames, a dynamic and visual display of the progress of multiple construction processes is achieved. The entire process realizes intelligent separation of changing signals, accurate calculation of engineering quantities, and realistic restoration of progress visualization. From statistical modeling and data measurement to visualization output, a full-link anti-interference mechanism is formed, greatly improving the accuracy of monitoring the progress of multiple construction processes in road engineering and the credibility of the visualization display.

[0171] Example 2, as Figure 4 As shown, based on the same inventive concept as the road construction progress visualization method based on UAV images provided in Embodiment 1, this embodiment of the invention also provides a road construction progress visualization system based on UAV images, including:

[0172] The aerial survey preprocessing module 11 is used to acquire raw images using a multispectral RTK UAV, obtain a corrected image set and identify potential semantic candidate blocks through image correction and road spatial stability block division;

[0173] The illumination decoupling module 12 is used to decompose global and local illumination interference through a multi-scale Retinex model based on a multispectral RTK UAV and a corrected image set, and decouple the pure material spectral reflectance component map and the real gradient component map to obtain the pure material spectral reflectance component map and the real gradient component map.

[0174] The feature fusion module 13 is used to retrieve the pure material spectral reflectance component map and the real gradient component map, construct an enhanced material spectral feature map set through block material spectral clustering attenuation weighted analysis, and fuse the neighborhood spectral gradient to construct a joint dataset that eliminates illumination anomalies.

[0175] Image registration module 14 is used to construct anchor point time series based on joint dataset and in combination with potential semantic candidate blocks. After Euclidean distance screening and cross-temporal spatial registration benchmark reconstruction, it obtains a time-stable semantic anchor point set and fits affine transformation to iteratively solve the exclusive affine transformation matrix of each period image.

[0176] The progress visualization module 15 is used to perform binary statistical feature separation modeling of construction changes and lighting disturbances based on a dedicated affine transformation matrix, obtain independent feature vectors, and use the effective pixel area as the statistical benchmark for engineering quantity to complete the visualization progress generation of complete time-series animation frame materials.

[0177] Furthermore, the aerial survey preprocessing module 11 is also used to: divide the aerial survey time sequence batches by road mileage, deploy encrypted flight routes on a multispectral RTK UAV, and record the optical imaging parameter set at each shooting time;

[0178] The original image is read and combined with the optical imaging parameter set. The radiance is solved pixel by pixel using the multispectral radiometric calibration extension algorithm. The radial distortion model is used to correct the pixel coordinate offset to obtain the corrected image set.

[0179] Each image in the corrected image set is divided into a road construction foreground area and an original landform background area using a fixed threshold coarse segmentation algorithm. The road CAD vector centerline is used as a spatial constraint to divide each image in the corrected image set into continuous blocks and align them with the centerline station number boundaries.

[0180] The average radiance of multiple periods is extracted block by block, and the temporal fluctuation difference is calculated to filter blocks whose fluctuations are lower than the preset amplitude threshold for multiple consecutive periods and are not within the road construction foreground area. These blocks are then marked as potential semantic candidate blocks.

[0181] Furthermore, the illumination decoupling module 12 is also used to: extract the solar elevation angle corresponding to each block, dynamically configure the Gaussian convolution kernel fusion weights, and use a multi-scale Retinex model to solve the global illumination component block by block, so as to obtain the intermediate radiation component map and the global illumination time series matrix stripped of global illumination disturbance.

[0182] Read the intermediate radiation component map corresponding to the current block to extract the average material reflection reference value of the potential semantic candidate block corresponding to the intermediate radiation component map;

[0183] By subtracting the average material reflection reference value from the intermediate radiance of each pixel within the block, the local illumination residual of a single pixel is obtained and a local illumination residual map is acquired.

[0184] By dividing the intermediate radiance by the pixel value at the same coordinate of the local illumination residual map through pixel-by-pixel division, the pure material spectral reflectance component map is obtained.

[0185] The pure material spectral reflectance component map is retrieved, and the true gradient component map is generated by performing pixel gradient calculation on the pure material spectral reflectance component map.

[0186] Furthermore, the feature fusion module 13 is also used to: read the mean value of global illumination components corresponding to all aerial survey periods for each block in the global illumination time series matrix, perform mean calculation on the mean value of global illumination components to obtain the baseline mean value, and subtract the baseline mean value from the mean value of global illumination components of the current single block to obtain the deviation difference value.

[0187] Retrieve the actual gradient component map corresponding to each block. If the pixel gradient is lower than the preset gradient threshold, it is determined to be an internal region of homogeneous material; otherwise, it is determined to be the boundary between different building materials.

[0188] Using the deviation difference as the input variable, a Gaussian attenuation weight function is constructed using a one-dimensional Gaussian distribution. For each pixel within the homogeneous material region, the Gaussian attenuation weight function is used for weighted correction to obtain the single-pixel corrected spectral value after offsetting the illumination shift and to construct a spectral feature map of the enhanced material.

[0189] Based on the real gradient component map, a standardized gradient mapping layer is obtained by linear normalization mapping of pixel gradient intervals.

[0190] Based on the standardized gradient mapping layer and the spectral feature atlas of the enhanced material, the mean value of the enhanced spectrum and the standardized gradient magnitude within the neighborhood are extracted and combined to form a binary joint vector, using a single pixel 3×3 neighborhood as the calculation unit.

[0191] The deviation difference is compared with the preset deviation threshold. If the deviation difference exceeds the preset deviation threshold, the corresponding spectral gradient joint vector is removed. After traversing all aerial survey time series batches, the effective binary joint vector is obtained.

[0192] The global illumination time series matrix, the deviation difference of each block, and the effective binary joint vector are combined to generate a joint dataset.

[0193] Furthermore, the image registration module 14 is also used to: based on the joint dataset, stitch together the single-pixel corrected spectral values ​​and pixel gradients of all pixels in each latent semantic candidate block to obtain anchor point binary vectors, and after extracting the anchor point binary vectors in the same latent semantic candidate block in a time-series batch, generate an anchor point time-series sequence.

[0194] Depend on Figure 3 It can be seen that the horizontal axis represents the index of the aerial survey period, and the vertical axis represents the Euclidean distance between the binary vectors of the anchor points of two adjacent periods;

[0195] Based on the anchor point time series, calculate the Euclidean distance between the binary vectors of two adjacent anchor points. Traverse all time series batches and record the potential semantic candidate blocks whose Euclidean distance is lower than the preset distance threshold for three consecutive periods as pre-stable anchor points.

[0196] Based on the mean Euclidean distance corresponding to each pre-stabilized anchor point, the anchor point temporal stability is graded and screened. Simultaneously, combined with the constraint of uniform distribution in the linear space of the road, the spatial distribution of anchor points is optimized to generate a set of temporally stable semantic anchor points.

[0197] The first phase of aerial survey imagery was selected as the reference imagery, and all other time-series aerial survey imagery were uniformly used as images to be registered.

[0198] Based on the temporally stable semantic anchor point set, the pixel coordinates of the corresponding stable semantic anchor points in the reference image and each image to be registered are extracted to construct coordinate pairs;

[0199] The affine transformation coordinate mapping equation is used to establish the pixel coordinate mapping relationship for the coordinate pairs. The coordinates of the stable semantic anchor points corresponding to the reference image are substituted into the affine transformation coordinate mapping equation to obtain the predicted coordinates. Combined with the pixel coordinates of the image to be registered in the corresponding coordinate pairs, the mean square of the deviation of all coordinate pairs is calculated and the square root is taken to obtain the root mean square error of the coordinates.

[0200] Using the root mean square error of coordinates as the loss function, several sets of coordinate pairs are substituted into the affine transformation coordinate mapping equation, and iterative calculations are performed until the loss function converges to a preset error threshold. The specific affine transformation matrix of the single image to be registered is then obtained.

[0201] Furthermore, the progress visualization module 15 is also used to: map the pure material spectral reflectance component map and the real gradient component map to the reference coordinate system based on the dedicated affine transformation matrix, solve the pure material spectral time difference and pixel gradient time difference of the same pixel in two adjacent periods pixel by pixel, so as to construct a binary difference vector, and simultaneously extract the illumination time difference at the same position to obtain a ternary statistical sample.

[0202] A binary Gaussian mixture model was built based on ternary statistical samples to classify two types of multivariate Gaussian distributions, including samples of illumination disturbance and samples of construction changes.

[0203] Using the difference in illumination time series as a prior constraint, the expectation-maximization algorithm is used to iteratively solve for the mean vector and covariance matrix corresponding to the two types of Gaussian distributions until the bivariate Gaussian mixture model converges iteratively. After convergence, the final mean vector and covariance matrix are obtained.

[0204] The optimal linear separating hyperplane is obtained based on the final mean vector and covariance matrix. All ternary statistical samples are substituted into the optimal linear separating hyperplane to complete the binary classification and generate a sample set with class labels.

[0205] Based on the sample set with attached category labels, extract the independent feature vectors corresponding to the two types of multivariate Gaussian distributions.

[0206] Extract the independent feature vectors corresponding to the construction change samples, and use the effective pixel area within the pixel region that conforms to the feature interval of the corresponding process as the statistical benchmark for the engineering quantity.

[0207] The pixel region corresponding to the process is traversed in each time sequence batch. The indicator function is used to determine whether a single pixel matches the feature interval of the corresponding process. If the match is successful, it is included in the effective construction area. If the match fails, it is not counted. The effective pixel areas of all collection periods before the current time sequence batch are accumulated in sequence and the area integration calculation is completed to obtain the cumulative engineering volume of the corresponding process up to the current time sequence batch.

[0208] Based on the cumulative workload, the current completion percentage of a single process is obtained, and the progress mapping function is obtained by weighted summation of multiple processes;

[0209] Based on the progress mapping function, using the pure material spectral reflectance component map and the real gradient component map as the visualization base, an independent rendering channel and layered progress color block are constructed for each process.

[0210] Independent rendering channels, layered progress color blocks, geographic anchor points, and station numbers are loaded in batches to generate visual animation frames, resulting in complete time-series animation frame materials to visualize the dynamic changes in the construction progress of multiple road construction processes throughout the entire time sequence.

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

[0212] This specification and accompanying drawings are merely illustrative examples of this application and are intended to cover any and all modifications, variations, combinations, or equivalents within the scope of this application. Clearly, those skilled in the art can make various alterations and modifications to this application without departing from its scope. Therefore, if such modifications and modifications fall within the scope of this application and its equivalents, this application intends to include such modifications and modifications.

Claims

1. A method for visualizing road construction progress based on UAV imagery, characterized in that, include: The original images were acquired using a multispectral RTK UAV. The corrected image set was obtained and potential semantic candidate blocks were identified through image correction and road spatial stability block division. Based on a multispectral RTK UAV and a corrected image set, global and local illumination interferences are decomposed hierarchically using a multi-scale Retinex model, and the pure material spectral reflectance component map and the true gradient component map are decoupled. By retrieving the pure material spectral reflectance component map and the true gradient component map, and constructing an enhanced material spectral feature map set through block material spectral clustering attenuation weighted analysis, and fusing the neighborhood spectral gradient to construct a joint dataset that eliminates illumination anomalies; Based on the joint dataset and the anchor point time series constructed by combining with the latent semantic candidate blocks, after Euclidean distance screening and cross-temporal spatial registration benchmark reconstruction, a time-stable semantic anchor point set is obtained, and the exclusive affine transformation matrix of each period image is solved by fitting affine transformation iteratively. Based on a dedicated affine transformation matrix, a binary statistical feature separation model of construction changes and lighting disturbances is performed to obtain independent feature vectors. The effective pixel area is used as the statistical benchmark for engineering quantity to complete the visualization progress generation of complete time-series animation frame materials.

2. The method for visualizing road construction progress based on UAV images according to claim 1, characterized in that, Original images were acquired using a multispectral RTK drone. Through image correction and road spatial stabilization segmentation, a corrected image set was obtained, and latent semantic candidate blocks were identified, including: The aerial survey was divided into batches based on road mileage, and a multispectral RTK drone was deployed to set up a dense flight path to record the optical imaging parameter set at each shooting moment; The original image is read and combined with the optical imaging parameter set. The radiance is solved pixel by pixel using the multispectral radiometric calibration extension algorithm. The radial distortion model is used to correct the pixel coordinate offset to obtain the corrected image set. Each image in the corrected image set is divided into a road construction foreground area and an original landform background area using a fixed threshold coarse segmentation algorithm. The road CAD vector centerline is used as a spatial constraint to divide each image in the corrected image set into continuous blocks and align them with the centerline station number boundaries. The average radiance of multiple periods is extracted block by block, and the temporal fluctuation difference is calculated to filter blocks whose fluctuations are lower than the preset amplitude threshold for multiple consecutive periods and are not within the road construction foreground area. These blocks are then marked as potential semantic candidate blocks.

3. The method for visualizing road construction progress based on UAV images according to claim 2, characterized in that, Based on a multispectral RTK UAV and a corrected image set, global and local illumination disturbances are decomposed hierarchically using a multi-scale Retinex model, resulting in the decoupled pure material spectral reflectance component map and the true gradient component map, including: The solar elevation angle corresponding to each block is extracted, and the Gaussian convolution kernel is dynamically configured to fuse the weights. The multi-scale Retinex model is used to solve the global illumination component block by block, and the intermediate radiation component map and global illumination temporal matrix after removing the global illumination perturbation are obtained. Read the intermediate radiation component map corresponding to the current block to extract the average material reflection reference value of the potential semantic candidate block corresponding to the intermediate radiation component map; By subtracting the average material reflection reference value from the intermediate radiance of each pixel within the block, the local illumination residual of a single pixel is obtained and a local illumination residual map is acquired. By dividing the intermediate radiance by the pixel value at the same coordinate of the local illumination residual map through pixel-by-pixel division, the pure material spectral reflectance component map is obtained. The pure material spectral reflectance component map is retrieved, and the true gradient component map is generated by performing pixel gradient calculation on the pure material spectral reflectance component map.

4. The method for visualizing road construction progress based on UAV images according to claim 3, characterized in that, The process of constructing a spectral feature atlas of reinforced materials includes: Read the mean global illumination component of each block in the global illumination time series matrix for all aerial survey periods, and perform mean calculation on the mean global illumination component to obtain the baseline mean. Subtract the baseline mean from the mean global illumination component of the current single block to obtain the deviation difference. Retrieve the actual gradient component map corresponding to each block. If the pixel gradient is lower than the preset gradient threshold, it is determined to be an internal region of homogeneous material; otherwise, it is determined to be the boundary between different building materials. Using the deviation difference as the input variable, a Gaussian attenuation weight function is constructed using a one-dimensional Gaussian distribution. For each pixel within the homogeneous material region, the Gaussian attenuation weight function is used for weighted correction to obtain the single-pixel corrected spectral value after offsetting the illumination shift and to construct a spectral feature map of the enhanced material.

5. The method for visualizing road construction progress based on UAV images according to claim 4, characterized in that, A joint dataset for removing illumination anomalies is constructed by fusing neighborhood spectral gradients, including: Based on the real gradient component map, a standardized gradient mapping layer is obtained by linear normalization mapping of pixel gradient intervals. Based on the standardized gradient mapping layer and the spectral feature atlas of the enhanced material, the mean value of the enhanced spectrum and the standardized gradient magnitude within the neighborhood are extracted and combined to form a binary joint vector, using a single pixel 3×3 neighborhood as the calculation unit. The deviation difference is compared with the preset deviation threshold. If the deviation difference exceeds the preset deviation threshold, the corresponding spectral gradient joint vector is removed. After traversing all aerial survey time series batches, the effective binary joint vector is obtained. The global illumination time series matrix, the deviation difference of each block, and the effective binary joint vector are combined to generate a joint dataset.

6. The method for visualizing road construction progress based on UAV images according to claim 5, characterized in that, Based on the joint dataset and combined with latent semantic candidate blocks to construct anchor point time series, after Euclidean distance filtering and cross-temporal spatial registration benchmark reconstruction, a time-stable semantic anchor point set is obtained, including: Based on the joint dataset, the single-pixel corrected spectral values ​​and pixel gradients of all pixels in each latent semantic candidate block are concatenated to obtain anchor binary vectors. After extracting anchor binary vectors in the same latent semantic candidate block in a time-series batch, an anchor time-series sequence is generated. Based on the anchor point time series, calculate the Euclidean distance between the binary vectors of two adjacent anchor points. Traverse all time series batches and record the potential semantic candidate blocks whose Euclidean distance is lower than the preset distance threshold for three consecutive periods as pre-stable anchor points. Based on the mean Euclidean distance corresponding to each pre-stabilized anchor point, the anchor point temporal stability is graded and screened. Simultaneously, combined with the constraint of uniform distribution in the linear space of the road, the spatial distribution of anchor points is optimized to generate a set of temporally stable semantic anchor points.

7. The method for visualizing road construction progress based on UAV images according to claim 6, characterized in that, The iterative solution of the affine transformation matrix for each period of the image is obtained by fitting the affine transformation, including: The first phase of aerial survey imagery was selected as the reference imagery, and all other time-series aerial survey imagery were uniformly used as images to be registered. Based on the temporally stable semantic anchor point set, the pixel coordinates of the corresponding stable semantic anchor points in the reference image and each image to be registered are extracted to construct coordinate pairs; The affine transformation coordinate mapping equation is used to establish the pixel coordinate mapping relationship for the coordinate pairs. The coordinates of the stable semantic anchor points corresponding to the reference image are substituted into the affine transformation coordinate mapping equation to obtain the predicted coordinates. Combined with the pixel coordinates of the image to be registered in the corresponding coordinate pairs, the mean square of the deviation of all coordinate pairs is calculated and the square root is taken to obtain the root mean square error of the coordinates. Using the root mean square error of coordinates as the loss function, several sets of coordinate pairs are substituted into the affine transformation coordinate mapping equation, and iterative calculations are performed until the loss function converges to a preset error threshold. The specific affine transformation matrix of the single image to be registered is then obtained.

8. The method for visualizing road construction progress based on UAV images according to claim 7, characterized in that, Based on a dedicated affine transformation matrix, a binary statistical feature separation model of construction changes and illumination disturbances is performed to obtain independent feature vectors, including: Based on the dedicated affine transformation matrix, the pure material spectral reflectance component map and the true gradient component map are mapped to the reference coordinate system. The pure material spectral time difference and pixel gradient time difference of the same pixel in two adjacent periods are solved pixel by pixel to construct a binary difference vector. Simultaneously, the illumination time difference at the same position is extracted to obtain a ternary statistical sample. A binary Gaussian mixture model was built based on ternary statistical samples to classify two types of multivariate Gaussian distributions, including samples of illumination disturbance and samples of construction changes. Using the difference in illumination time series as a prior constraint, the expectation-maximization algorithm is used to iteratively solve for the mean vector and covariance matrix corresponding to the two types of Gaussian distributions until the bivariate Gaussian mixture model converges iteratively. After convergence, the final mean vector and covariance matrix are obtained. The optimal linear separating hyperplane is obtained based on the final mean vector and covariance matrix. All ternary statistical samples are substituted into the optimal linear separating hyperplane to complete the binary classification and generate a sample set with class labels. Based on the sample set with attached category labels, extract the independent feature vectors corresponding to the two types of multivariate Gaussian distributions.

9. The method for visualizing road construction progress based on UAV images according to claim 8, characterized in that, Using the effective pixel area as the benchmark for engineering quantity statistics, the visualization progress of complete time-series animation frame materials is generated, including: Extract the independent feature vectors corresponding to the construction change samples, and use the effective pixel area within the pixel region that conforms to the feature interval of the corresponding process as the statistical benchmark for the engineering quantity. The pixel region corresponding to the process is traversed in each time sequence batch. The indicator function is used to determine whether a single pixel matches the feature interval of the corresponding process. If the match is successful, it is included in the effective construction area. If the match fails, it is not counted. The effective pixel areas of all collection periods before the current time sequence batch are accumulated in sequence and the area integration calculation is completed to obtain the cumulative engineering volume of the corresponding process up to the current time sequence batch. Based on the cumulative workload, the current completion percentage of a single process is obtained, and the progress mapping function is obtained by weighted summation of multiple processes; Based on the progress mapping function, using the pure material spectral reflectance component map and the real gradient component map as the visualization base, an independent rendering channel and layered progress color block are constructed for each process. Independent rendering channels, layered progress color blocks, geographic anchor points, and station numbers are loaded in batches to generate visual animation frames, resulting in complete time-series animation frame materials to visualize the dynamic changes in the construction progress of multiple road construction processes throughout the entire time sequence.

10. A road construction progress visualization system based on UAV imagery, used to implement the road construction progress visualization method based on UAV imagery as described in any one of claims 1 to 9, characterized in that, include: The aerial survey preprocessing module is used to acquire raw images using a multispectral RTK UAV, obtain a corrected image set and identify potential semantic candidate blocks through image correction and road spatial stability block division; The illumination decoupling module is used to decompose global and local illumination interference through a multi-scale Retinex model based on a multispectral RTK UAV and a corrected image set, and decouple the pure material spectral reflectance component map and the real gradient component map. The feature fusion module is used to retrieve the pure material spectral reflectance component map and the true gradient component map. Through the attenuation weighted analysis of the block material spectral clusters, it constructs an enhanced material spectral feature map set and fuses the neighborhood spectral gradient to construct a joint dataset that eliminates illumination anomalies. The image registration module is used to construct a temporal sequence of anchor points based on the joint dataset and in combination with latent semantic candidate blocks. After Euclidean distance filtering and cross-temporal spatial registration benchmark reconstruction, a temporally stable set of semantic anchor points is obtained, and the affine transformation matrix of each period image is solved iteratively by fitting affine transformation. The progress visualization module is used to perform binary statistical feature separation modeling of construction changes and lighting disturbances based on a dedicated affine transformation matrix, obtain independent feature vectors, and use the effective pixel area as the statistical benchmark for engineering quantity to complete the visualization progress generation of complete time-series animation frame materials.

Images