Land environment shadow area optical feature recovery method and system

Abstract

The invention discloses a land environment shadow region optical feature recovery method and system, and the method comprises the steps: synchronously collecting a visible light image, a near infrared spectrum image and three-dimensional point cloud data of a target land region through a multi-mode sensor, and generating multi-source fusion data with a space-time alignment attribute; based on the multi-source fusion data, constructing an illumination transmission physical model of the scene, and separating a direct illumination attenuation component and an indirect illumination diffuse reflection component of a shadow region; performing compensation reconstruction on the direct illumination attenuation component by using a generative adversarial network, and generating an initial recovery image after shadow elimination; and outputting a final recovered image with global illumination consistency and local detail authenticity by adopting a cross-scale attention fusion mechanism. By using the embodiment of the invention, the real spectral attribute and rich spatial texture details of the ground feature can be maintained, and the visual consistency of the image and the reliability of subsequent quantitative application are improved.

Description

Technical Field

[0001] This invention belongs to the field of remote sensing technology, and in particular to a method and system for restoring optical features of shadowed areas in a land environment. Background Technology

[0002] In the field of land remote sensing, shadowed areas are prevalent in optical images due to factors such as topographic relief, cloud cover, and building obstruction. These shadows severely distort the true spectral reflectance characteristics of ground features, leading to a significant decrease in the accuracy of subsequent image-based applications such as land feature classification, change detection, and quantitative inversion. Traditional shadow processing methods often rely on a single temporal phase or a single data source, making it difficult to accurately model the radiative transfer processes of complex surfaces under non-uniform illumination. This often results in color distortion, loss of detail, or inconsistencies in illumination with surrounding non-shadowed areas after shadow removal. Although some studies have attempted to improve the situation by combining multispectral or lidar data, limitations remain in the synergistic recovery of spectral fidelity and spatial details in shadowed areas. The lack of a physical separation and targeted compensation mechanism for direct illumination and indirect illumination from multiple reflections further hinders the realization of high-precision optical feature recovery. Summary of the Invention

[0003] The purpose of this invention is to provide a method and system for restoring optical features of shadowed areas in terrestrial environments, in order to overcome the shortcomings of the prior art, maintain the true spectral properties and rich spatial texture details of ground features, and improve the visual consistency of images and the reliability of subsequent quantitative applications.

[0004] One embodiment of this application provides a method for restoring optical features of shadowed areas in a land environment, the method comprising: Visible light images, near-infrared spectral images, and 3D point cloud data of the target land area are acquired simultaneously by multimodal sensors to generate multi-source fusion data with spatiotemporal alignment attributes; Based on the multi-source fusion data, a physical model of light transmission in the scene is constructed. By solving the reflectivity equation under non-uniform illumination field, the direct light attenuation component and the indirect light diffuse reflection component in the shadow area are separated. Generative adversarial networks are used to compensate and reconstruct the direct illumination attenuation component, and spectral consistency correction is performed based on the indirect illumination diffuse reflection component to generate an initial restored image after shadow removal. A cross-scale attention fusion mechanism is adopted to adaptively weight and fuse the initial restored image with the original highlight region features, and output the final restored image with global illumination consistency and local detail authenticity.

[0005] Another embodiment of this application provides a system for restoring optical features of shadowed areas in a terrestrial environment, the system comprising: The acquisition module is used to simultaneously acquire visible light images, near-infrared spectral images, and three-dimensional point cloud data of the target land area through multimodal sensors, and generate multi-source fusion data with spatiotemporal alignment attributes; The construction module is used to construct a physical model of light transmission in the scene based on the multi-source fusion data, and to separate the direct light attenuation component and the indirect light diffuse reflection component in the shadow area by solving the reflectivity equation under non-uniform light field. The reconstruction module is used to compensate and reconstruct the direct illumination attenuation component using a generative adversarial network, and to perform spectral consistency correction based on the indirect illumination diffuse reflection component, thereby generating an initial restored image after shadow removal. The output module is used to adaptively weight and fuse the initial restored image with the original highlight region features using a cross-scale attention fusion mechanism, and output the final restored image with global illumination consistency and local detail authenticity.

[0006] Another embodiment of this application provides an electronic device including a memory and a processor, wherein the memory stores a computer program and the processor is configured to run the computer program to perform the method described in any of the preceding claims.

[0007] Compared with existing technologies, the optical feature restoration method for shadowed areas of land environment provided by the present invention can maintain the true spectral properties and rich spatial texture details of ground objects, improve the visual consistency of images and the reliability of subsequent quantitative applications. Attached Figure Description

[0008] Figure 1 A hardware structure block diagram of a computer terminal for a method of restoring optical features of shadowed areas in a terrestrial environment, provided in an embodiment of the present invention; Figure 2 A flowchart illustrating a method for restoring optical features of shadowed areas in a terrestrial environment, provided in an embodiment of the present invention; Figure 3 This is a schematic diagram of the structure of a system for restoring optical features of shadowed areas in a terrestrial environment, provided in an embodiment of the present invention. Detailed Implementation

[0009] The embodiments described below with reference to the accompanying drawings are exemplary and are only used to explain the present invention, and should not be construed as limiting the present invention.

[0010] This invention first provides a method for restoring optical features of shadowed areas in a terrestrial environment. This method can be applied to electronic devices, such as computer terminals, specifically ordinary computers.

[0011] The following detailed explanation uses a computer terminal as an example. Figure 1This is a hardware structure block diagram of a computer terminal for a method of restoring optical features of shadowed areas in a terrestrial environment, provided in an embodiment of the present invention. Figure 1 As shown, the computer device includes a processor, memory, and network interface connected via a system bus, wherein the memory may include non-volatile storage media and internal memory.

[0012] See Figure 2 The present invention provides a method for restoring optical features of shadowed areas in a land environment, which may include the following steps: S201 uses multimodal sensors to simultaneously acquire visible light images, near-infrared spectral images, and three-dimensional point cloud data of the target land area, generating multi-source fusion data with spatiotemporal alignment attributes; Specifically, visible light cameras, near-infrared spectrometers, and lidar sensors can be integrated on the same observation platform. Hardware triggering circuits ensure that all sensors start acquiring data at the same time, generating raw synchronous multimodal data streams. The core of this step is to achieve the physical integration and synchronous acquisition of multimodal sensors, and to ensure the spatiotemporal consistency of data through hardware-level triggering. The specific implementation method is as follows: The selection of the observation platform must be suitable for the terrestrial environment observation scenario. Stable mobile platforms or fixed observation frames should be given priority. The platform must have a sensor mounting reference surface with a flatness error ≤0.05mm, ensuring that the parallelism deviation of the sensor optical axis is ≤0.1° to avoid spatial acquisition deviations caused by installation tilt. When integrating sensors, the visible light camera, near-infrared spectrometer, and lidar sensor should be fixed on the reference surface, with their optical axes aligned and the center-to-center distance controlled within the range of 10-20cm to reduce parallax effects. At the same time, sufficient space should be reserved for the wiring of sensor data lines and trigger lines to avoid signal distortion caused by line interference.

[0013] The hardware trigger circuit design is the core of synchronous data acquisition. It employs a master-slave trigger architecture, using a high-precision clock module as the trigger signal source. The clock frequency is set to 100MHz (to ensure trigger signal accuracy, with a period error ≤1ns). The main trigger signal, generated by the trigger circuit, is transmitted to the trigger interfaces of the three sensors via three synchronous trigger lines. The trigger signal uses a pulse signal with the following parameters: pulse amplitude 5V (adapting to the sensor trigger voltage range of 3.3-5V), pulse width 10μs (ensuring stable sensor recognition of the trigger signal), and trigger frequency 10Hz (i.e., acquiring 10 frames of data per second, balancing acquisition efficiency and data volume). To ensure trigger synchronization, shielded cables are used, with consistent cable lengths (error ≤5cm) to avoid signal transmission delay differences. The synchronization error is controlled within ≤1μs, ensuring all sensors start exposure or data acquisition at the same time.

[0014] The process of generating the original synchronous multimodal data stream is as follows: When the hardware trigger circuit outputs a trigger pulse, the visible light camera immediately starts exposure, with the exposure time set to 1 / 100s (to adapt to daytime land observation scenarios and avoid overexposure or underexposure), generating RGB image data with a resolution of 1920×1080; the near-infrared spectrometer simultaneously starts spectral acquisition, with the acquisition wavelength range set to 700-1100nm and the spectral resolution of 5nm, generating spectral image data containing 200 spectral channels, with the image resolution consistent with that of the visible light camera (1920×1080); the lidar sensor simultaneously starts laser emission and reception, with a laser emission frequency of 100kHz, a scanning angle range of 0-360°, and a vertical resolution of 0.1°, generating three-dimensional point cloud data, each point cloud containing three-dimensional coordinates (x, y, z) and reflection intensity information. All three raw data streams are accompanied by a collection timestamp with a precision of 1μs and a format of "YYYY-MM-DDHH:MM:SS.ssssss". They are ultimately aggregated into a raw synchronous multimodal data stream and stored in the platform's local storage module. The data transmission rate is ≥100Mbps to ensure no data loss.

[0015] The original synchronous multimodal data stream is timestamped and its spatial coordinate system is unified. The pre-calibrated sensor extrinsic parameter matrix is ​​used to transform the data of each mode to the same spatial coordinate system, generating multimodal data with unified spatiotemporal reference. The core of this step is to eliminate the temporal deviation and spatial coordinate system differences of the original data and establish a unified spatiotemporal reference. The specific implementation method is as follows: Timestamp correction addresses potential minute time offsets in the original data (caused by sensor response delays) using linear interpolation. First, the timestamp sequences of the three sensors are extracted from the original data stream. Let the timestamp sequence of the visible light camera be T1=[t1_1,t1_2,...,t1_n], the near-infrared spectrometer be T2=[t2_1,t2_2,...,t2_n], and the lidar be T3=[t3_1,t3_2,...,t3_n]. Using the lidar timestamp as a reference (where lidar timestamp stability is optimal), the time deviations between T1, T2, and T3 are calculated as Δt1=T1-T3 and Δt2=T2-T3. If the deviation Δt ≤ 1μs, it is considered to be in good synchronization and no correction is required; if the deviation 1μs < Δt ≤ 10μs, the timestamp of the data frame is adjusted by linear interpolation. For example, the timestamp of a visible light image frame t1_k = 1699999999.000005s corresponds to the timestamp of the lidar t3_k = 1699999999.000003s, with a deviation Δt1 = 2μs. The corrected timestamp t1_k' = t3_k + (t1_k - t3_k) × k / n (k is the current frame number, and n is the total number of frames), ensuring that the timestamps of all modal data are completely consistent after correction.

[0016] The core of Spatial Coordinate System 1 is to achieve coordinate transformation through pre-calibrated extrinsic parameter matrices. First, the lidar coordinate system is selected as the unified spatial coordinate system (right-handed coordinate system, x-axis pointing in the observation direction, y-axis horizontal and perpendicular to the x-axis, z-axis vertically upward). The pre-calibration process is completed using a checkerboard calibration board, resulting in the extrinsic parameter matrix M1 (3×4 matrix, including rotation matrix R1 (3×3) and translation vector T1 (3×1)) of the visible light camera relative to the lidar, and the extrinsic parameter matrix M2 (3×4 matrix, including rotation matrix R2 and translation vector T2) of the near-infrared spectrometer relative to the lidar. The physical meaning of the extrinsic parameter matrices is to describe the attitude relationship between the sensor coordinate system and the unified coordinate system. For example, the rotation matrix R1 is used to correct the sensor's installation angle deviation, and the translation vector T1 is used to correct the sensor's spatial position offset. The pre-calibration accuracy requires a rotation error ≤ 0.1° and a translation error ≤ 1cm.

[0017] The coordinate transformation process is as follows: For each pixel (u, v) in the visible light image, its corresponding 3D coordinates (xc, yc, zc) in the camera coordinate system are converted to world coordinates through the camera intrinsic parameter matrix (pre-calibrated, including parameters such as focal length f and principal point coordinates (u0, v0)). Then, it is converted to unified coordinate system coordinates (x, y, z) = R1 × (xc, yc, zc) + T1 through the extrinsic parameter matrix M1. The pixel coordinates of the near-infrared spectral image undergo the same process, using the intrinsic and extrinsic parameter matrices M2 to convert them to unified coordinate system coordinates. The original point cloud data from the lidar is already in its own coordinate system and is directly used as the original data for the unified coordinate system. For example, for a visible light pixel (u=960, v=540), after intrinsic parameter transformation, the camera coordinates are (0, 0, 10)m, the extrinsic parameter matrix R1 is the identity matrix (without rotational bias), and T1 is (0.1, 0, 0)m. After conversion, the unified coordinate system coordinates are (10.1, 0, 0)m, achieving unification with the lidar point cloud coordinates. The final result is multimodal data with a unified spatiotemporal reference, in which each pixel or point cloud is associated with a unified timestamp and spatial coordinates.

[0018] Based on multimodal data with unified spatiotemporal reference, a feature point matching algorithm is used to accurately register 3D point clouds with visible light images and near-infrared spectral images, realizing a one-to-one correspondence between pixels and point clouds, and generating registered multi-source data; The core of this step is to eliminate spatial parallax between sensors through feature point matching, and establish a precise correlation between pixels and point clouds. The specific implementation method is as follows: The feature point matching algorithm employs a scale-invariant feature transform algorithm, which possesses scale and rotation invariance and can effectively handle differences in viewpoints between sensors. The core algorithm process includes feature point extraction, feature descriptor generation, feature matching, and mismatch removal. In the feature point extraction stage, Gaussian difference pyramids are constructed for both visible light and near-infrared spectral images. The pyramid has 6 layers, and the Gaussian difference scale factor is 1.6. Candidate feature points are obtained through extremum detection, followed by edge response removal and sub-pixel localization to obtain a stable set of feature points. For example, approximately 2000 feature points can be extracted from a 1920×1080 visible light image, each containing position (u,v), scale, and orientation information.

[0019] In the feature descriptor generation stage, a 16×16 neighborhood window is selected centered on each feature point. This window is then divided into 4×4 sub-windows. Gradient histograms in eight directions are calculated for each sub-window, ultimately generating a 128-dimensional feature descriptor. The Euclidean distance of the descriptor is used to measure the similarity between two feature points. For 3D point clouds, feature points are extracted by calculating the point cloud's normal vectors and curvature (e.g., using the SIFT-3D algorithm), generating a 128-dimensional descriptor with the same dimension as the image feature points, thus achieving matchability between the point cloud and image feature points.

[0020] During feature matching, a nearest neighbor matching strategy is employed, with a distance threshold of 0.6. This means that a valid match is considered complete when the ratio of the nearest neighbor descriptor distance to the second nearest neighbor distance for a given image feature point is less than 0.6. For example, feature point A in a visible light image has a Euclidean distance of 0.3 to feature point B in the point cloud and a distance of 0.5 to feature point C in the point cloud. Since the ratio 0.6 ≤ 0.6, A and B are considered a valid match. To further improve matching accuracy, a random sampling consensus algorithm is used to eliminate false matches. The algorithm iterates 1000 times, with an inlier threshold of 2cm (meaning a spatial distance of less than 2cm between matching points is considered an inlier). The final matching accuracy is ≥95%.

[0021] Establishing a one-to-one correspondence between pixels and point clouds: By effectively matching feature point pairs, a mapping model between image pixels and point clouds is constructed (using homography matrix or perspective n-point algorithm). This model can map any pixel (u,v) in the image to the corresponding point (x,y,z) in the 3D point cloud. For example, by training a perspective n-point model with 1000 effective matching point pairs, for any pixel (1000,600) in the image, the corresponding point cloud coordinates (15.2,3.1,0.8)m can be calculated through the model, achieving a one-to-one correspondence between all image pixels and point clouds. The registration process between near-infrared spectral images and point clouds is consistent with that of visible light images, ultimately generating registered multi-source data. Each visible light pixel and near-infrared pixel in the data is associated with unique 3D point cloud coordinates and reflection intensity information.

[0022] Multi-dimensional feature fusion is performed on the registered multi-source data to extract RGB channel features of visible light images, reflectance features of near-infrared spectral images, and geometric elevation features of three-dimensional point clouds, generating multi-source fused data with spatiotemporal alignment attributes.

[0023] The core of this step is to integrate the core features of multimodal data to form a unified feature representation. The specific implementation method is as follows: Multidimensional feature extraction requires ensuring the consistency and fusionability of the dimensions of each feature. RGB channel feature extraction from visible light images: The grayscale values ​​of the R, G, and B channels of each pixel are normalized using the formula f_rgb=(I-I_min) / (I_max-I_min), where I is the original grayscale value (0-255), I_min=0, I_max=255, and the normalized feature value range is [0,1]. A 3D RGB channel feature vector is generated for each pixel. For example, a pixel with original R=200, G=150, B=100 will have a normalized feature vector of [0.784, 0.588, 0.392], which reflects the color information of the target area.

[0024] Reflectance feature extraction from near-infrared spectral images: The raw data acquired by the near-infrared spectrometer are spectral intensity values, which need to be converted into reflectance values. The conversion formula is R = I_sample / I_standard, where I_sample is the spectral intensity of the target area, and I_standard is the spectral intensity of the standard white board (pre-acquired and stored). The reflectance R ranges from [0,1]. Each pixel corresponds to the reflectance values ​​of 200 spectral channels, generating a 200-dimensional reflectance feature vector. This feature reflects the spectral characteristics of the material in the target area. For example, if the intensity value of a pixel in the 700nm channel is 1200, and the corresponding channel intensity of the standard white board is 2000, then the reflectance of this channel is 0.6. The reflectance values ​​of the 200 channels constitute the reflectance feature vector of this pixel.

[0025] Geometric elevation feature extraction of 3D point clouds: Using the z-axis coordinate of a unified coordinate system as the elevation feature, the normal vector (x, y, z components) and curvature of each point cloud are calculated simultaneously to generate a 4D geometric elevation feature vector (z, nx, ny, nz), where z is the elevation value (in meters), nx, ny, and nz are the three components of the normal vector (range [-1, 1]), and curvature is the degree of bending of the point cloud surface (range [0, 1]). For example, if a point cloud has coordinates (15.2, 3.1, 0.8) m, a normal vector of (0.02, 0.01, 0.9997), and a curvature of 0.03, then the geometric elevation feature vector is [0.8, 0.02, 0.01, 0.9997, 0.03]. This feature reflects the spatial geometry of the target area.

[0026] Multidimensional feature fusion employs a feature concatenation method, sequentially concatenating the 3D RGB features, 200-dimensional near-infrared reflectance features, and 5-dimensional geometric elevation features corresponding to each pixel to generate a 208-dimensional multi-source fusion feature vector. During the fusion process, it is crucial to ensure that all features are associated with the same timestamp and spatial coordinates to achieve spatiotemporal alignment. For example, the RGB features of a pixel [0.784, 0.588, 0.392], near-infrared reflectance features (200 values ​​ranging from 0.5 to 0.8), and geometric elevation features [0.8, 0.02, 0.01, 0.9997, 0.03], after concatenation, form a 208-dimensional vector. Each vector is labeled with its corresponding timestamp (e.g., 1699999999.000003s) and spatial coordinates (15.2, 3.1, 0.8) m. The final multi-source fusion data is stored in the form of a set of feature vectors, with the data format being "timestamp-spatial coordinates-208-dimensional fusion feature", ensuring that the information of each modality can be accurately associated in subsequent processing.

[0027] S202, Based on the multi-source fusion data, construct a physical model of light transmission in the scene, and separate the direct light attenuation component and the indirect light diffuse reflection component in the shadow area by solving the reflectivity equation under non-uniform light field. Specifically, the three-dimensional geometric structure, surface normal vector distribution, and material reflectivity characteristics of the scene can be extracted from multi-source fusion data. Combined with solar azimuth angle, elevation angle, and atmospheric transmittance parameters, a physical model of light transmission under non-uniform illumination field can be constructed. The core of this step is to extract the core physical attributes of the scene from multi-source fusion data, integrate external lighting and atmospheric parameters, and establish a physical model that can accurately describe the transmission law of light in the scene. The specific implementation method is as follows: The extraction of the scene's 3D geometry is based on the 3D point cloud coordinate information from multi-source fusion data, and is accomplished using the Poisson reconstruction algorithm. First, the 3D point cloud (x, y, z) coordinates associated with each pixel are separated from the multi-source fusion data to form the original point cloud dataset. The point cloud density is set to 100 points / square centimeter (to ensure complete depiction of scene details, such as ground texture and low-lying vegetation outlines). The Poisson reconstruction algorithm constructs an indicator function for the point cloud and uses the Poisson equation to solve for a smooth 3D surface. The algorithm parameters are set as follows: reconstruction depth level 8 (corresponding to a reconstruction resolution of 0.1 meters, capable of distinguishing scene bumps at the 10-centimeter level), sampling density threshold 0.05 (to remove sparse noise points), and surface smoothness factor 0.8 (balancing detail preservation and surface smoothness). In the example, the target land area includes urban roads, roadside trees and low shrubs. After reconstruction, the original point cloud can clearly show the flat surface of the road (z-coordinate fluctuation ≤ 0.05 meters), the cylindrical outline of the tree trunk (diameter about 0.3 meters) and the irregular volume shape of the shrubs. The generated three-dimensional mesh model has about 500,000 vertices, which completely restores the spatial geometry of the scene.

[0028] The extraction of surface normal vector distribution is achieved through point cloud neighborhood analysis. For the reconstructed 3D mesh model, the k nearest neighbors of each vertex (k=20, balancing computational efficiency and normal vector accuracy) are taken to construct the covariance matrix of the neighborhood points. After eigenvalue decomposition of the covariance matrix, the eigenvector corresponding to the smallest eigenvalue is the surface normal vector of that vertex. The direction of the normal vector is corrected by the overall scene coordinate system (z-axis vertically upward) to ensure that the normal vector points outward from the scene. In the example, the z-coordinates of the neighborhood points of a vertex on the road surface are all within the range of 0.8-0.85 meters, and the vector corresponding to the smallest eigenvalue of the covariance matrix is ​​(0.02, 0.01, 0.9997), which is corrected to (0.02, 0.01, 0.9997), indicating that the surface of this area is close to horizontal; the neighborhood points of a vertex on the tree trunk are distributed in a ring, and the normal vector is (0.98, 0.03, 0.18), pointing radially outward from the tree trunk, accurately reflecting the cylindrical surface characteristics of the tree trunk. The final generated surface normal vector distribution is attached to the 3D mesh model in the form of a vector field, with each vertex associated with a 3D normal vector component (range [-1, 1]) of (nx, ny, nz).

[0029] The extraction of material reflectance characteristics fuses visible light and near-infrared spectral features. First, the 3D RGB normalized features and 200-dimensional near-infrared reflectance features of each pixel are separated from the multi-source fused data. A spectral angle matching algorithm is then used to classify materials and calibrate their reflectance. The algorithm calculates the angle between the target pixel's spectrum and a standard material spectral library; a smaller angle indicates higher similarity, with a threshold of 5° (identifying them as the same material). The standard material spectral library includes common terrestrial materials (concrete, asphalt, vegetation, and soil). For example, the near-infrared reflectance of concrete is stable at 0.3-0.4 in the 800-1100nm band, while the reflectance of vegetation abruptly jumps to above 0.6 in the 800nm ​​band (vegetation red-edge effect). In the example, the near-infrared spectrum of a certain pixel has a reflectance of 0.62 in the 800nm ​​band, and the angle between it and the standard spectrum of vegetation is 3.2°. It is determined to be vegetation material, and its material reflectance characteristics are determined as follows: visible light RGB reflectance (0.15, 0.4, 0.2) and near-infrared average reflectance 0.65. This characteristic directly characterizes the material's ability to reflect light in different wavelength bands.

[0030] Acquisition and quantification of solar azimuth, altitude, and atmospheric transmittance parameters: The solar azimuth is defined as the angle between the projection of sunlight onto the horizontal plane and true north (range 0-360°), and the solar altitude is the angle between sunlight and the horizontal plane (range 0-90°). These are calculated using the GPS coordinates (latitude 30°N, longitude 120°E) and time (Beijing time 10:00). In this example, the azimuth is 35° and the altitude is 42°. Atmospheric transmittance characterizes the degree of light attenuation after passing through the atmosphere. It is calculated based on the atmospheric window bands (940nm, 1064nm) of the near-infrared spectrum, using the formula τ=exp(-k×d), where k is the atmospheric extinction coefficient (taken as 0.05 / km, based on the standard atmospheric model), and d is the observation distance (2km). The calculated τ=exp(-0.05×2)=0.9048, meaning that after 2km of atmospheric transmission, the intensity of the sunlight is 90.48% retained.

[0031] The physical model of light transmission under non-uniform illumination field is based on the theory of radiative transmission. The core formula is: L(x,ω)=ρ(x)×[E_dir(x,ω_s)×τ×cosθ(x)+E_ind(x)] / π+L_sky(x,ω). The parameters are as follows: L(x,ω) is the radiance of pixel x in the observation direction ω (unit: W / (m²・sr)); ρ(x) is the material reflectivity of pixel x (extracted from the material reflectivity characteristics); E_dir(x,ω_s) is the direct irradiance in the solar direction ω_s (taken as 1000W / m² under standard clear weather conditions); τ is the atmospheric transmittance; θ(x) is the angle between the sunlight and the surface normal vector of pixel x (calculated from the solar azimuth angle, altitude angle and surface normal vector, cosθ(x)=nx×sinh×cosaz+ny×sinh×sinaz+nz×cosh, where h is the solar altitude angle and az is the azimuth angle); E_ind(x) is the indirect irradiance (the light intensity reflected from the surrounding environment, initially taken as 100W / m²); L_sky(x,ω) is the sky radiance (taken as 50W / (m²・sr)). In the example, for the road pixel x, ρ(x) = 0.3 (asphalt material), θ(x) = 30° (cosθ = 0.866), and τ = 0.9048. Substituting these values ​​into the formula, we get L(x,ω) = 0.3 × [1000 × 0.9048 × 0.866 + 100] / π + 50 ≈ 0.3 × [783.5 + 100] / 3.1416 + 50 ≈ 0.3 × 883.5 / 3.1416 + 50 ≈ 84.5 + 50 = 134.5 W / (m²・sr), accurately depicting the illumination radiation characteristics of this pixel. The model adapts to the illumination differences in different areas through non-uniformity correction terms (the variation of E_dir and E_ind with scene position), such as the attenuation of E_dir in areas under tree shade and the enhancement of E_ind in areas obscured by buildings.

[0032] Based on the physical model of light transmission, a reflectivity equation for the shadow region is established. This equation describes the attenuation process of direct light after it is blocked and the contribution of diffuse reflection formed by indirect light through multiple reflections. The core of this step is to define the range of the shadow area, establish a reflectivity equation specifically, and quantify the contributions of direct light attenuation and indirect light diffuse reflection. The specific implementation method is as follows: The definition of shadow areas is determined by the occlusion relationship between the lighting direction and the 3D geometry. This is based on the 3D geometric model constructed in step one and the solar direction vector (determined by an azimuth angle of 35° and an elevation angle of 42°, with the vector being (sin42°×cos35°, sin42°×sin35°, cos42°) ≈ (0.544, 0.368, 0.743)). For each pixel x in the scene, the solar direction vector is traced in reverse. If other geometric objects (such as tree trunks or buildings) occlude the tracing path (occlusion distance < distance from the observation point to the sun), then x is determined to be a shadow area. In the example, the projection of the tree trunks (x=15m, y=3m, z=3m) along the solar direction vector covers the road area (x=16-18m, y=3.2-3.4m, z=0.8m). This road area is determined to be a shadow area, with a shadow edge transition width of approximately 0.2 meters (consistent with real-world lighting diffraction effects).

[0033] The reflectance equation is established based on the physical model of light transmission. Addressing the core characteristic of shadow regions—the attenuation of direct light due to shading—a direct light attenuation coefficient α(x) is introduced (α(x)∈[0,1], α(x)=1 in non-shaded areas, α(x)<1 in shadowed areas, representing the proportion of direct light retention). Simultaneously, considering the contribution of diffuse reflection from indirect light, an indirect light diffuse reflection coefficient β(x) is introduced (β(x)∈[0,1], representing the proportion of diffuse reflection from the surrounding environment to this area). The final reflectance equation for the shadow region is: ρ_obs(x)=ρ_true(x)×[α(x)×E_dir×τ×cosθ(x)+β(x)×∫ΩE_env(ω')×ρ_env(ω')×cosθ'(x)dω'] / E_ref. The parameters are as follows: ρ_obs(x) is the observed reflectance of pixel x in the shadow region (converted from the grayscale value of the visible light image, the conversion formula is ρ_obs(x)=I(x) / (I_max), where I(x) is the grayscale value and I_max=255); ρ_true(x) is the true material reflectance of pixel x (extracted in step one, unaffected by illumination); ∫ΩE_env(ω')×ρ_env(ω')×cosθ'(x)dω' is the indirect illumination integral term, Ω is the hemispherical space (all possible reflected light directions), E_env(ω') is the ambient light intensity, ρ_env(ω') is the ambient material reflectance, θ'(x) is the angle between the ambient reflected light and the surface normal vector of pixel x; E_ref is the reference illumination intensity (taken as 1000W / m², uniform dimension).

[0034] Quantification of direct light attenuation process: Obstacles (such as tree trunks) block direct sunlight, causing the direct light intensity in the shadow area to decrease proportionally to the projected area of ​​the obstacle. The initial value of α(x) is determined by the projection coverage of the obstacle at pixel x, and the formula is α(x) = 1 - S_occ(x) / S_pix, where S_occ(x) is the projected area of ​​the obstacle on the plane where pixel x is located, and S_pix is ​​the actual ground area corresponding to the pixel (in the example, the pixel resolution is 0.1m × 0.1m, S_pix = 0.01m², and the shadow pixel S_occ(x) = 0.008m², then α(x) = 1 - 0.008 / 0.01 = 0.2, that is, only 20% of the direct light is retained). During the attenuation process, the supplementary effect of atmospheric scattering also needs to be considered, and α(x) needs to be corrected. The correction formula is α'(x)=α(x)+(1-α(x))×τ_scatter, where τ_scatter is the atmospheric scattering coefficient (taken as 0.15, characterizing the proportion of the scattered light supplementing the shadow area). After correction, α'(x)=0.2+0.8×0.15=0.32.

[0035] Quantification of diffuse reflection contribution from indirect lighting: Indirect lighting originates from multiple reflections from other non-shadowed areas in the scene. For example, the non-shadowed ground around a shadowed area reflects direct lighting into the shadowed area, forming a diffuse reflection contribution. The integral term is calculated using the Monte Carlo sampling method, randomly sampling 100 lighting directions ω' within the hemispherical space Ω. The E_env(ω') for each direction is calculated from the ambient pixel reflectivity ρ_env(ω') and the direct lighting intensity (E_env(ω') = ρ_env(ω') × E_dir × τ × cosθ_env, where θ_env is the angle between the sunlight and the ambient pixel surface normal vector). In the example, the surrounding environment of a pixel x in the shaded area is mainly a concrete road surface (ρ_env=0.35). The average E_env(ω') of 100 directions is 280W / m², and the average cosθ'(x) is 0.6. Therefore, the integral term result is 280×0.6=168W / m². β(x) is taken as 0.8 (characterizing that 80% of the integral term illumination can reach pixel x), so the indirect illumination contribution term is 0.8×168=134.4W / m².

[0036] Example of the physical meaning of the reflectance equation: The observed reflectance of a certain shadow pixel x is ρ_obs(x)=0.12 (grayscale value I(x)=31), the real material reflectance is ρ_true(x)=0.3 (asphalt pavement), E_dir=1000W / m², τ=0.9048, cosθ(x)=0.7, α(x)=0.32, the indirect illumination integral term is 168W / m², β(x)=0.8, and E_ref=1000W / m². Substituting into the equation for verification: Right side = 0.3×[0.32×1000×0.9048×0.7+0.8×168] / 1000=0.3×[202.7+134.4] / 1000=0.3×337.1 / 1000=0.101. The deviation from the observed reflectance of 0.12 is caused by noise, which will be corrected through subsequent algorithm optimization.

[0037] Using visible light images and near-infrared spectral images from multi-source fusion data as real observations, the least squares optimization algorithm is used to solve the reflectance equation to obtain the direct illumination attenuation coefficient and indirect illumination diffuse reflectance coefficient of the shadow region. The core of this step is to establish a correlation between image observations and the reflectance equation, and to solve for key coefficients through least squares optimization to adapt the equation to the real scene. The specific implementation method is as follows: Acquisition and preprocessing of real observations: The visible light image and near-infrared spectral image of the shadow area are separated from the multi-source fusion data. The RGB gray values ​​of the visible light image are converted into the observed reflectance ρ_obs_rgb(x) (red, green and blue bands). The conversion formula is ρ_obs_rgb(x)=I_rgb(x) / (E_ref×G), where G is the camera gain (taken as 0.001W / (m²・sr・DN), obtained by camera calibration). In the example, the R channel gray value of a certain shadow pixel is I_r=25, then ρ_obs_r(x)=25 / (1000×0.001)=0.025. Three characteristic bands (750nm, 850nm, and 950nm, covering the characteristic spectral regions of vegetation and soil) were selected from the near-infrared spectral image and converted into the observed reflectance ρ_obs_nir(x). The conversion formula is the same as that for visible light. In the example, the gray value I_nir of the 850nm band is 40, so ρ_obs_nir(x) = 40 / (1000×0.001) = 0.04. Finally, each shadow pixel obtained 6 observation values ​​(RGB + 3 near-infrared bands), forming the observation vector Y(x) = [ρ_obs_r(x), ρ_obs_g(x), ρ_obs_b(x), ρ_obs_nir1(x), ρ_obs_nir2(x), ρ_obs_nir3(x)].

[0038] The objective function of the least squares optimization algorithm is constructed as follows: Based on the reflectance equation, the residual relationship between the observed values ​​and the model prediction values ​​is established, and the optimization objective is to minimize the sum of squared residuals of all shadow pixels. The model prediction value ρ_pred(x) is calculated by the reflectance equation, i.e., ρ_pred(x) = ρ_true(x) × [α × E_dir × τ × cosθ(x) + β × E_ind(x)] / E_ref, where α (direct illumination attenuation coefficient) and β (indirect illumination diffuse reflection coefficient) are variables to be optimized (global variables; α and β have consistent values ​​within the same shadow area to ensure the continuity of illumination distribution), and ρ_true(x), E_dir, τ, cosθ(x), and E_ind(x) are all known quantities (extracted or calculated in steps one and two). The objective function is: min(Σ||Y(x_i)-ρ_pred(x_i;α,β)||²), where x_i is the i-th pixel in the shadow region (1000 pixels are selected to construct an optimization sample set, balancing computational efficiency and accuracy), and ||・|| is the L2 norm.

[0039] The constraints for the optimization process are set as follows: α∈[0,1] (the attenuation coefficient cannot be negative or greater than 1, and its physical meaning is the light retention ratio), β∈[0,1] (the diffuse reflection coefficient is similar), to avoid the optimization results having values ​​with unreasonable physical meaning. The least squares optimization is solved using the Gauss-Newton iterative method. The iterative steps are as follows: 1. Initialize the optimization variables α0=0.3 and β0=0.7 (based on empirical values); 2. Calculate the predicted value ρ_pred(x_i;α0,β0) and residual r_i=Y(x_i)-ρ_pred(x_i;α0,β0) for each sample; 3. Construct the Jacobian matrix J (a 6000×2 matrix, 6 observations × 1000 samples, each row corresponding to the partial derivative of an observation with respect to α and β); 4. Solve the increment equation J^TJΔx=-J^Tr (Δx=[Δα,Δβ] is the variable increment); 5. Update the variables α1=α0+Δα and β1=β0+Δβ; 6. Calculate the residual sum of squares. If the residual sum of squares is less than the threshold (1e-6) or the number of iterations reaches 50, stop the iteration; otherwise, repeat steps 2-5.

[0040] Example optimization process: Initially, α0=0.3 and β0=0.7, and the initial residual sum of squares is calculated to be 0.08; in the first iteration, Δα=0.02 and Δβ=0.03 are obtained. After updating, α1=0.32 and β1=0.73, and the residual sum of squares is reduced to 0.03; after 15 iterations, the residual sum of squares is reduced to 0.0008 (less than the threshold 1e-6), and the iteration is stopped. Finally, the optimized direct illumination attenuation coefficient α=0.35 and indirect illumination diffuse reflection coefficient β=0.76 are obtained. Optimization result verification: Substituting α=0.35 and β=0.76 into the reflectance equation, the predicted observation vector for a certain sample pixel is [0.024,0.031,0.028,0.039,0.042,0.037]. The average deviation from the actual observation vector [0.025,0.032,0.029,0.040,0.043,0.038] is 0.001, which is less than 5%, meeting the accuracy requirements. This result indicates that only 35% of the direct illumination in this shaded area is retained, while the diffuse reflection contribution of indirect illumination accounts for 76%, consistent with the illumination distribution pattern of tree-shaded shadows.

[0041] Based on the attenuation coefficient and diffuse reflection coefficient obtained from the solution, the image data of the shadow area is decomposed into a direct illumination attenuation component and an indirect illumination diffuse reflection component, generating two independent illumination component images.

[0042] The core of this step is to reverse-decompose the illumination contribution of the shadow region image based on the optimized coefficients, thereby achieving the separation and visualization of the two illumination components. The specific implementation method is as follows: Preprocessing of image data in the shadow area: The visible light image (resolution 1920×1080) of the shadow area is selected as the decomposition object. The image grayscale values ​​are converted into radiance L(x) (unit: W / (m²・sr)). The conversion formula is L(x)=I(x)×G+B, where B is the camera dark current offset (taken as 0.1W / (m²・sr) to eliminate the influence of camera noise). In the example, the grayscale value I(x) of a certain shadow pixel is 30, then L(x)=30×0.001+0.1=0.13W / (m²・sr). This radiance includes the contribution of direct illumination attenuation and the contribution of indirect illumination diffuse reflection, which is the basis for subsequent decomposition.

[0043] The decomposition formulas for the direct illumination attenuation component and the indirect illumination diffuse reflection component are derived from the radiance form of the reflectivity equation: Radiance L(x) = L_dir(x) + L_ind(x), where L_dir(x) is the direct illumination attenuation component (the radiance of direct illumination after attenuation due to shading), and L_ind(x) is the indirect illumination diffuse reflection component (the radiance of ambient diffuse illumination). Combining the optimization coefficients α and β from step three, the decomposition formulas are: L_dir(x) = ρ_true(x) × α × E_dir × τ × cosθ(x) / π; L_ind(x) = ρ_true(x) × β × E_ind(x) / π. The meanings of the parameters in the formula are consistent with those in the previous text: ρ_true(x) is the true reflectivity of the material, α is the direct illumination attenuation coefficient, β is the indirect illumination diffuse reflection coefficient, E_dir is the direct illumination intensity, τ is the atmospheric transmittance, cosθ(x) is the cosine value of the incident angle, E_ind(x) is the indirect illumination intensity, and π is pi (a unit of uniform radiance).

[0044] Example of decomposition process: For a certain shadow pixel x, ρ_true(x)=0.3 (asphalt material), α=0.35, E_dir=1000W / m², τ=0.9048, cosθ(x)=0.7, β=0.76, E_ind(x)=168W / m². Substituting into the formula, we can calculate: L_dir(x)=0.3×0.35×1000×0.9048×0.7 / 3.1416≈0.3×0.35×633.36 / 3.1416≈0.3×70.37 / 3.1416≈6.7W / (m²・sr); L_ind(x)=0.3×0.76×168 / 3.1416≈0.3×127.68 / 3.1416≈12.1W / (m²・sr). The total radiance of the pixel is L(x) = 6.7 + 12.1 = 18.8 W / (m²·sr), which deviates from the preprocessed L(x) = 18.7 W / (m²·sr) by only 0.1 W / (m²·sr), and the decomposition accuracy meets the requirements.

[0045] Generation of two independent illumination component images: For each pixel in the shadow area, calculate L_dir(x) and L_ind(x) sequentially, convert the radiance value into image grayscale value (conversion formula I_dir(x)=L_dir(x)×255 / L_max, L_max=50W / (m²・sr), which is the maximum radiance of the scene), and generate the direct illumination attenuation component image and the indirect illumination diffuse reflection component image. The characteristics of the direct illumination attenuation component image are: low gray values ​​in the shadow core region (small L_dir(x), e.g., I_dir(x) = 30-50), gradually increasing gray values ​​at the shadow edges (due to reduced occlusion, increased α, and increased L_dir(x), e.g., I_dir(x) = 80-120), and the highest gray values ​​in the non-shadow region (α = 1, I_dir(x) = 150-200). The characteristics of the indirect illumination diffuse reflection component image are: uniform gray value distribution (small fluctuations in E_ind(x) within the shadow region), gray value range of 50-80, with a slight increase only near buildings or dense vegetation (enhanced ambient diffuse reflection). The generated two component images have the same resolution as the original visible light image (1920×1080), with pixel positions corresponding one-to-one, which can intuitively present the distribution differences of the two illumination components in the shadow region, providing accurate component data for subsequent illumination compensation.

[0046] S203, using a generative adversarial network to compensate and reconstruct the direct illumination attenuation component, and performing spectral consistency correction based on the indirect illumination diffuse reflection component, to generate an initial restored image after shadow removal; Specifically, a generative adversarial network model can be constructed, in which the generator takes the direct illumination attenuation component as input, the discriminator takes the non-shaded area image as reference, and learns the illumination compensation mapping relationship through adversarial training to obtain a trained generative adversarial network. The core of this step is to build a generative adversarial network architecture adapted to the illumination compensation task. Through adversarial game between the generator and the discriminator, the mapping relationship from the direct illumination attenuation component to the normal illumination component is learned. The specific implementation method is as follows: The overall architecture of the generative adversarial network model adopts a dual-module collaborative structure of generator and discriminator. The generator is responsible for illumination compensation and reconstruction, while the discriminator is responsible for real and fake image discrimination. The generator uses the U-Net architecture, which has a symmetrical encoder-decoder structure and skip connections, effectively preserving image detail information and adapting to the fine reconstruction requirements of illumination components. The encoder consists of six convolutional layers, each with the following parameters: kernel size 3×3, stride 2, padding 1, and LeakyReLU activation function (negative slope 0.2, enhancing gradient propagation). This downsamples and extracts features from the input direct illumination attenuation component, progressively downsampling from a 64×64 resolution image to a 1×1 resolution, generating a 512-dimensional deep feature vector. The decoder consists of six deconvolutional layers, each with the same parameters: kernel size 3×3, stride 2, padding 1, and ReLU activation function. This upsamples the deep feature vector to a 64×64 resolution and simultaneously fuses shallow detail features (such as shadow edge textures) from the corresponding encoder layers via skip connections. The final output is a compensated reconstructed direct illumination component with the same size as the input.

[0047] The discriminator uses the PatchGAN architecture, which focuses on distinguishing between real and fake images in local regions and is better suited to the requirements of local illumination consistency in illumination compensation tasks. The discriminator contains four convolutional layers with the following parameters: kernel size 4×4, stride 2, padding 1, and LeakyReLU activation function (negative slope 0.2). The last layer uses the Sigmoid activation function to output a 30×30 patch-level probability map (each patch corresponds to a 70×70 local region of the input image), with probability values ​​∈ [0,1], used to determine whether the input image is a "fake image" generated by the generator or a "real image" of a real non-shaded region.

[0048] The construction of the training dataset must ensure the consistency of data distribution. The dataset contains 10,000 pairs of samples of the target land area. Each pair of samples consists of a "direct illumination attenuation component image" and a "corresponding non-shaded image of the area". The sample collection covers different land scenes (urban roads, farmland, woodland) and different lighting conditions (sunny days, cloudy days) to ensure the model's generalization ability. The data preprocessing process includes: scaling the image resolution to a uniform 64×64, normalizing the gray values ​​to the range of [0,1] (the normalization formula is I_norm=(I-I_min) / (I_max-I_min), where I is the original gray value, I_min=0, I_max=255), and randomly horizontally flipping and rotating (rotation angle ±10°) to enhance data diversity.

[0049] The adversarial training process employs an alternating training strategy. Training parameters are set as follows: batch size 16, initial learning rate 0.0002 (using the Adam optimizer with momentum parameters β1=0.5 and β2=0.999 to improve training stability), and a total of 50,000 iterations. The training loss function uses a combination of adversarial loss and L1 reconstruction loss. The adversarial loss is calculated using the cross-entropy between the generator and discriminator, with the formula L_adv=-E[log(D(x_true))]-E[log(1-D(G(x_input)))], where D is the discriminator, G is the generator, x_true is the real non-shadow image, and x_input is the direct illumination attenuation component. The L1 reconstruction loss constrains the pixel-level error between the generated and real images, with the formula L1=||G(x_input)-x_true||1 / N, where N is the total number of pixels. The combined loss weight is L_total=0.1×L_adv+L1. During training, the model parameters are saved every 1000 iterations. Training stops when the combined loss value fluctuates less than 1e-4 for 1000 consecutive iterations. In the example, L_total=0.023 after 38000 iterations, reaching the convergence state and obtaining the trained generative adversarial network.

[0050] The direct illumination attenuation component is input into the trained generative adversarial network to generate the compensated and reconstructed direct illumination component, which restores the illumination intensity that should be in the shadow area. The core of this step is to use the trained generator model to compensate for the direct illumination attenuation component in the shadow area, restoring its normal illumination intensity in the unobstructed state. The specific implementation method is as follows: The preprocessing of the input data must be consistent with that of the training phase to ensure that the model input format matches. First, from the direct illumination attenuation component image generated in step four, image blocks of the shadow area (resolution 64×64) are cropped, and the gray value of each pixel is normalized. In the example, the original gray value of a certain shadow pixel is 45 (corresponding to weak illumination intensity, which is consistent with the characteristics of the shadow area), and after normalization, it is (45-0) / 255≈0.176. At the same time, in order to eliminate the influence of image noise on the compensation effect, Gaussian filtering is used to preprocess the input image. The filter kernel size is 3×3, the standard deviation σ=0.5, and the filtering formula is G(x,y)=1 / (2πσ²)×exp(-(x²+y²) / (2σ²)). After filtering the pixel (5,5), the gray value is corrected from 45 to 44, and the noise fluctuation is reduced by 30%.

[0051] The compensation and reconstruction process is implemented by calling a pre-trained generator model. The pre-processed direct illumination attenuation component image patch is input into the generator. The model extracts illumination attenuation features (such as grayscale distribution and edge gradients) from the input image through the encoder, and then fuses these features through upsampling and skip connections in the decoder, finally outputting the compensated and reconstructed direct illumination component image patch. The generator's output needs to be denormalized to restore it to the original grayscale value range. The denormalization formula is I_restore=I_gen×(I_max-I_min)+I_min, where I_gen is the normalized grayscale value output by the model. In the example, the normalized grayscale value output by the model is 0.72. After denormalization, I_restore=0.72×255≈183, corresponding to a significant increase in illumination intensity, approaching the normal illumination level of non-shaded areas (average grayscale value of non-shaded areas is 185).

[0052] The compensation effect was verified through quantitative analysis of light intensity. Light intensity and grayscale value have a linear relationship (light intensity I_light = k × I_restore, where k is the conversion coefficient, calibrated to k = 0.01 W / (m²・sr・DN)). In the example, before compensation, the light intensity in the shadow area was I_light_pre = 45 × 0.01 = 0.45 W / (m²・sr), and after compensation, I_light_post = 183 × 0.01 = 1.83 W / (m²・sr). The deviation from the light intensity of the adjacent non-shadow area (1.85 W / (m²・sr)) was only 1.08%, achieving the expected restoration of light intensity in the shadow area. Simultaneously, the detail preservation effect was verified using an edge detection algorithm (Canny algorithm, threshold set to 50-150). The compensated image clearly retains details such as ground cracks and vegetation textures, with the difference in edge gradient values ​​between the image and the non-shadow area ≤ 5%, showing no detail blurring or distortion.

[0053] Using the diffuse reflection component of indirect illumination as the spectral reference, the spectral difference between the direct illumination component and the adjacent non-shaded area after compensation and reconstruction is calculated, and a spectral correction algorithm is used to eliminate spectral distortion in the compensation process. The core of this step is to correct the spectral distortion that may occur during the generative adversarial network compensation process based on the stable spectral characteristics of the diffuse reflection component of indirect illumination, thus ensuring the spectral consistency after illumination compensation. The specific implementation method is as follows: The spectral reference is selected based on the diffuse reflection component of indirect illumination because it is dominated by ambient diffuse reflection and is not affected by direct sunlight. Its spectral characteristics are highly consistent with the true spectrum of the material, providing a stable reference value. From the diffuse reflection component image generated in step four, the spectral data of the non-shaded area within a 50-pixel range adjacent to the shaded area are extracted as the spectral reference set for that shaded area. The spectral reference set includes reflectance data for the three visible light channels (R: 650nm, G: 550nm, B: 450nm) and three near-infrared characteristic channels (750nm, 850nm, 950nm). The reflectance calculation formula is ρ = I_restore / (I_standard), where I_standard is the grayscale value of the standard white board (pre-acquired at 255). In the example, the reflectance data of the spectral reference set are R: 0.32, G: 0.45, B: 0.28, 750nm: 0.52, 850nm: 0.63, and 950nm: 0.58, which correspond to the typical spectral characteristics of vegetation materials (the reflectance in the 850nm band is significantly higher than that in the visible light band).

[0054] The spectral difference is calculated using a spectral angle matching algorithm. This algorithm quantifies spectral similarity by calculating the angle between two spectral vectors; the smaller the angle, the better the spectral consistency. First, the reflectance of the six feature bands corresponding to the compensated and reconstructed direct illumination component image is extracted as a spectral vector V_gen=[ρ_gen_R,ρ_gen_G,ρ_gen_B,ρ_gen_750,ρ_gen_850,ρ_gen_950]. The average reflectance of the spectral reference set is extracted as a reference vector V_ref=[ρ_ref_R,ρ_ref_G,ρ_ref_B,ρ_ref_750,ρ_ref_850,ρ_ref_950]. The formula for calculating the spectral angle θ is cosθ=(V_gen・V_ref) / (||V_gen||×||V_ref||), where ・ is the dot product and ||・|| is the L2 norm. In the example, V_gen=[0.25,0.38,0.22,0.48,0.55,0.53], V_ref=[0.32,0.45,0.28,0.52,0.63,0.58], and we can calculate cosθ≈0.978 and θ≈12.8°. We set the spectral angle threshold to 5° (exceeding the threshold indicates spectral distortion), and spectral correction is required.

[0055] The spectral correction algorithm employs a polynomial fitting correction method, selecting a second-order polynomial as the fitting function to establish a mapping relationship between the compensated reflectance and the reference reflectance. The fitting function is ρ_corr = a × ρ_gen² + b × ρ_gen + c, where ρ_corr is the corrected reflectance, and a, b, and c are the fitting coefficients. The coefficients are solved by least-squares fitting of the spectral reference set and the compensated spectral data. The fitting objective is to minimize the sum of squared residuals between the corrected reflectance and the reference reflectance.

[0056] The spectrally corrected direct illumination component is linearly superimposed with the indirect illumination diffuse reflection component to generate an initial restored image after shadow removal, which has eliminated shadows but retains the original texture details.

[0057] The core of this step is to restore the complete lighting information of the shadow area by linearly superimposing and fusing the two lighting components, generating an initial restored image with shadow removal and detail preservation. The specific implementation method is as follows: The weighting coefficients for linear superposition are determined based on the energy conservation principle of the physical model of light transmission. The two light components contribute differently to the total radiance of the image. The direct light component is the primary source of illumination, and its weighting coefficient is set to ω_dir = 0.7. The indirect diffuse light component is a supplementary source of illumination, and its weighting coefficient is set to ω_ind = 0.3, with ω_dir + ω_ind = 1 (to ensure energy conservation). The weighting coefficients can be verified by the proportion of scene illumination intensity. In non-shaded areas, the direct illumination intensity accounts for approximately 72%, and the indirect illumination accounts for approximately 28%. The error between the set weighting coefficients and the actual illumination distribution is ≤2%, ensuring the physical rationality of the superposition result.

[0058] The linear superposition formula is I_init(x,y)=ω_dir×I_corr_dir(x,y)+ω_ind×I_ind(x,y), where I_init(x,y) is the gray value of the initial restored image at pixel (x,y), I_corr_dir(x,y) is the gray value of the direct illumination component after spectral correction, and I_ind(x,y) is the gray value of the diffuse reflection component under indirect illumination. Before superposition, it is necessary to ensure that the pixel positions of the two component images are completely aligned (based on the spatiotemporal alignment in step two), and the gray value range is uniformly [0,255]. In the example, for pixel (200, 300), I_corr_dir=183 and I_ind=62. Substituting these values ​​into the formula, we get I_init=0.7×183+0.3×62=128.1+18.6=146.7, which is rounded to 147. This grayscale value is within the range of the average grayscale value (145-150) of the non-shadow area, thus eliminating the shadow effect.

[0059] The detail preservation verification during the overlay process is achieved through texture feature extraction. The gray-level co-occurrence matrix (GLCM) algorithm is used to extract texture feature parameters (contrast, correlation, energy, entropy) from the image, and the differences in texture parameters between the initial restored image and the original non-shaded area image are compared. The parameters of the GLCM are set as follows: distance d=1, angle θ=0°, 45°, 90°, 135°, and 16 gray levels are counted. In the example, the contrast of the original non-shaded area is 85, the correlation is 0.82, the energy is 0.15, and the entropy is 1.8; the contrast of the initial restored image is 83, the correlation is 0.81, the energy is 0.14, and the entropy is 1.78. The difference in each parameter is ≤2.5%, indicating that the overlay image completely preserves the original texture details, such as the texture of gravel on the ground and the veins of leaves in vegetation.

[0060] The post-processing of the initial restored image includes contrast enhancement and noise removal. An adaptive histogram equalization algorithm (CLAHE algorithm, clipLimit=2.0, tileGridSize=8×8) is used to enhance local contrast and improve the visual effect of the image. A median filtering algorithm (3×3 kernel size) is used to remove minor noise that may be generated during the overlay process. The median filtering formula is I_final(x,y)=median{I_init(xi,yj)|i,j∈[-1,1]}. In the example, the grayscale values ​​of pixels surrounding pixel (200,300) are 145, 147, 146, 148, 147, 149, 146, 148, and 147, with a median of 147. After filtering, the grayscale value remains unchanged, removing noise without blurring details. The final generated initial restored image has the same resolution as the original visible light image (1920×1080). Shadow areas are completely eliminated, illumination distribution is uniform, global illumination difference is ≤3%, and local texture details are clear, meeting the basic requirements for subsequent fusion optimization.

[0061] S204. A cross-scale attention fusion mechanism is adopted to adaptively weight and fuse the initial restored image with the original highlight region features, and output the final restored image with global illumination consistency and local detail authenticity.

[0062] Specifically, multi-scale features of the highlight region can be extracted from the original visible light image, and feature pyramids of the corresponding scale can be extracted from the initial restored image to generate a highlight region feature set and an initial restored feature pyramid. The core of this step is to simultaneously extract multi-scale features from two types of images, establish a feature alignment basis, and provide a suitable feature carrier for subsequent cross-scale fusion. The specific implementation method is as follows: Defining the highlight region is a prerequisite for feature extraction. A dual judgment method combining brightness threshold and local contrast is adopted. First, the gray values ​​of the original visible light image are normalized to the range of [0, 255]. A brightness threshold T = 220 is set (based on the statistics of clear daylight in land scenes, areas above this value are potential highlight regions), and a set of pixels with brightness greater than T is initially screened out. Then, the local contrast C of this pixel set is calculated. The local contrast formula is C = (I_max - I_min) / I_mean, where I_max and I_min are the maximum and minimum gray values ​​in a 3×3 neighborhood, respectively, and I_mean is the neighborhood mean. A contrast threshold C_th = 0.2 is set. When C ≥ C_th, it is determined to be a real highlight region (excluding bright spot noise).

[0063] Multi-scale feature extraction of highlight regions is achieved using a backbone network based on convolutional neural networks. A network architecture with strong feature representation capabilities is selected, comprising five convolutional stages. Each stage downsamples through convolutional and pooling layers, generating feature maps at four scales (corresponding to 1 / 1, 1 / 2, 1 / 4, and 1 / 8 of the original image resolution, forming a multi-scale hierarchy). The convolutional layer parameters are uniformly set as follows: kernel size 3×3, stride 1, padding 1, and ReLU activation function to ensure non-linear representation of feature extraction. Max pooling is used in the pooling layers with a kernel size of 2×2 and a stride of 2, achieving downsampling while preserving key features. For highlight regions, the highlight mask is multiplied pixel-by-pixel with the feature maps from each stage, retaining only the feature information of the highlight region to generate multi-scale highlight features. In the example, the original image resolution is 1920×1080, and the specular feature map sizes at the four scales are 1920×1080×64, 960×540×128, 480×270×256, and 240×135×512 respectively (the number of channels doubles as the scale decreases, improving the ability to express deep features). The pixel value of each feature map corresponds to the texture, edge and other detailed features of the specular region, forming a specular region feature set.

[0064] The initial reconstruction feature pyramid extraction uses the same backbone network and parameter settings as the specular feature extraction to ensure scale alignment. The initial reconstruction image is input into the backbone network and processed through five convolutional stages, generating four feature maps at corresponding scales, with dimensions identical to the specular feature set (1920×1080×64 to 240×135×512), forming the initial reconstruction feature pyramid. The core function of the feature pyramid is to capture image information at different scales. Shallow features (1 / 1 and 1 / 2 scales) focus on local texture details (such as ground cracks and vegetation textures), while deep features (1 / 4 and 1 / 8 scales) focus on global illumination distribution and semantic information (such as region material types). In the example, the 1 / 2 scale feature map (960×540×128) of the initial reconstruction image clearly captures the texture edge features of the road after shadow removal, while the 1 / 8 scale feature map (240×135×512) characterizes the illumination distribution trend of the entire scene, providing multi-dimensional feature support for subsequent cross-scale fusion.

[0065] Design a cross-scale attention fusion module that calculates the similarity between highlight region features and initial restored features at different scales and generates a scale-aware attention weight matrix. The core of this step is to adaptively capture the correlation between two types of features at different scales through an attention mechanism, generate an accurate weight matrix, and achieve "on-demand fusion". The specific implementation method is as follows: The overall architecture of the cross-scale attention fusion module includes four sub-modules: feature alignment, similarity calculation, scale-aware weighting, and weight normalization. The feature alignment sub-module first unifies the dimensions of the highlight features and the initial restored features at different scales. For the k-th scale (k=1,2,3,4 corresponding to 4 scales), let the highlight feature be F_h^k∈R^(H_k×W_k×C_k), and the initial restored feature be F_r^k∈R^(H_k×W_k×C_k) (H_k and W_k are the height and width of the feature map, and C_k is the number of channels). The number of channels of the two types of features is uniformly mapped to C=256 through a 1×1 convolutional layer (reducing the computational cost while preserving the feature representation). The mapping formula is F_h^k'=Conv1×1(F_h^k) and F_r^k'=Conv1×1(F_r^k), where the number of kernels in the Conv1×1 convolution is 256, the stride is 1, and the padding is 0. In the example, the F_h^2 of the second scale (960×540×128) is convolved by 1×1 to generate F_h^2'∈R^(960×540×256), which has the same dimension as F_r^2' of the same scale, thus completing feature alignment.

[0066] The similarity calculation submodule uses cosine similarity to quantify the similarity between two types of aligned features. Cosine similarity can effectively measure the directional consistency of feature vectors, which meets the requirements for association determination of illumination features. For each pixel position (i,j) at each scale, the cosine similarity s_ij^k between F_h^k'(i,j) and F_r^k'(i,j) is calculated, with the formula s_ij^k=(F_h^k'(i,j)・F_r^k'(i,j)) / (||F_h^k'(i,j)||×||F_r^k'(i,j)||), where ・ is the dot product, ||・|| is the L2 norm, and s_ij^k∈[-1,1]. The closer the value is to 1, the more similar the features are, and the closer it is to -1, the greater the difference. In the example, the dot product of F_h^3'(i,j) and F_r^3'(i,j) of a pixel (i=100,j=200) at the third scale (480×270×256) is 125.6, and the L2 norms are 15.2 and 14.8 respectively. The calculated s_ij^3=125.6 / (15.2×14.8)≈125.6 / 224.96≈0.558 indicates that there is a moderate degree of similarity between the two types of features.

[0067] The scale-aware weighted submodule is used to balance the similarity weights at different scales. Considering that shallow features (large scale) focus on details and deep features (small scale) focus on the global picture, a scale weight coefficient λ_k is introduced. λ_k increases with the scale level k (λ_1=0.2 when k=1, λ_2=0.3 when k=2, λ_3=0.25 when k=3, λ_4=0.25 when k=4, and the sum is 1), achieving a balance between detailed and global features. The scale-aware similarity is s_ij^k'=λ_k×s_ij^k. Then, the similarity is mapped to the range [0,1] by the Sigmoid function to obtain the initial attention weight w_ij^k=1 / (1+exp(-s_ij^k')). The closer w_ij^k is to 1, the more highlight features need to be fused at that position; the closer it is to 0, the more the initial recovery features are preserved.

[0068] The weight normalization submodule performs row normalization on the initial weight matrix to ensure that the sum of the weights at each pixel position is 1, avoiding local weight imbalance. The normalization formula is w_ij^k'=w_ij^k / Σ(i,j)w_ij^k, which ultimately generates a scale-aware attention weight matrix W^k∈R^(H_k×W_k×1) for each scale. The weight matrices of the four scales form a weight set, which corresponds one-to-one with the scale of the feature pyramid, providing accurate weight guidance for subsequent weighted fusion.

[0069] The initial restored feature pyramid is weighted and fused using an attention weight matrix, and the detailed features of the highlight region are adaptively fused into the initial restored image to generate a fused feature map. The core of this step is to achieve adaptive fusion of two types of features based on attention weights. While preserving the illumination consistency of the initial restored image, it accurately integrates the detailed features of the highlight areas. The specific implementation method is as follows: The core logic of weighted fusion is "highlight feature guidance + initial feature dominance". The fusion formula is designed for the feature map at each scale k as: F_fuse^k(i,j)=w_ij^k'×F_h^k'(i,j)+(1-w_ij^k')×F_r^k'(i,j), where F_fuse^k is the fusion feature at scale k, and w_ij^k' is the normalized attention weight. This formula achieves an adaptive effect of "more fusion of similar regions and less fusion of different regions" through weight allocation. For regions with high similarity between highlight features and initial restored features (w_ij^k'→1), more highlight details are incorporated; for regions with large differences (w_ij^k'→0), the illumination consistency of the initial restored image is preserved first to avoid illumination distortion introduced by highlight features.

[0070] Example fusion process: For a pixel (i=500, j=800) at scale 1 (1920×1080×256), w_ij^1'=0.8 (high similarity, more highlights need to be fused), the mean eigenvector of F_h^1'(i,j) is 0.65 (corresponding to a region with rich highlight details), and the mean eigenvector of F_r^1'(i,j) is 0.42 (initial recovered texture features). After fusion, F_fuse^1(i,j)=0.8×0.65+0.2×0.42=0 0.52 + 0.084 = 0.604, which retains the basic texture of the initial restoration while incorporating the detailed features of the highlights; for another pixel (i=300, j=400), w_ij^1'=0.1 (low similarity, less specular fusion), F_h^1'(i,j) mean 0.7, F_r^1'(i,j) mean 0.38, after fusion = 0.1×0.7 + 0.9×0.38 = 0.07 + 0.342 = 0.412, mainly retaining the initial restored features and avoiding specular interference.

[0071] During the fusion process, it is necessary to ensure the spatial alignment of feature maps at different scales. The downsampling / upsampling stride is strictly controlled by the pooling and convolution parameters of the backbone network to ensure that the pixel positions of feature maps at different scales can be accurately mapped to the coordinates of the original image. For example, the fused feature map F_fuse^4 at the 4th scale (240×135×256) has each pixel corresponding to an 8×8 region of the original image, which can be accurately restored to the original resolution through subsequent upsampling.

[0072] Multi-scale fusion feature aggregation: The fusion feature maps F_fuse^1 to F_fuse^4 at four scales are aggregated across scales. Upsampling and stitching are used. First, the deep small-scale features (F_fuse^3, F_fuse^4) are upsampled to the same size as the shallow large-scale features (F_fuse^1, F_fuse^2) through bilinear interpolation. Then, the feature maps of the same size are stitched along the channel dimension to generate the aggregated fusion feature. In the example, F_fuse^4 (240×135×256) is upsampled twice by bilinear interpolation (scaled by a factor of 2 each time), resulting in a size of 960×540×256. This size is then concatenated with F_fuse^2 (960×540×256) to generate a 960×540×512 aggregated feature map. F_fuse^3 (480×270×256) is upsampled once to become 960×540×256. This aggregated feature map is then concatenated with the above aggregated feature map to generate a 960×540×768 feature map. Finally, this feature map is upsampled to 1920×1080×768 and concatenated with F_fuse^1 (1920×1080×256) to generate a final fused feature map of 1920×1080×1024. This feature map integrates the specular details and initial recovery features of four scales, taking into account both local details and global illumination consistency.

[0073] The initial verification of the fusion effect was achieved through feature entropy calculation. A higher entropy value indicates richer feature information, and the entropy value of the fused feature map needs to be higher than that of the initial restored feature map. In the example, the entropy value of the initial restored feature map was 1.85, and the entropy value of the fused feature map was 2.32, an improvement of 25.4%, indicating that the highlight detail features have been effectively integrated, and the feature information is richer.

[0074] The fused feature map is decoded and reconstructed, the image resolution is restored through deconvolution, and the local texture is enhanced by a detail enhancement algorithm. The final output is a restored image with global illumination consistency and local detail realism.

[0075] The core of this step is to decode the high-dimensional fusion features into a visual image. Through resolution restoration and detail enhancement, it balances global illumination uniformity and local texture clarity. The specific implementation method is as follows: Decoding and reconstruction are achieved using a deconvolutional network. This network forms a symmetrical structure with the backbone network used for feature extraction described earlier. It consists of five deconvolutional stages, corresponding to feature recovery at five scales. Each deconvolutional stage comprises a deconvolutional layer, a batch normalization layer, and a ReLU activation function. The deconvolutional layer parameters are set as follows: kernel size 3×3, stride 2, padding 1, ensuring that each deconvolution can double the feature map size (doubling the resolution) and halve the number of channels (reducing computational cost). For example, the input fused feature map is 240×135×1024 (deep aggregated features). After 5 deconvolutions, the size becomes 480×270×512, 960×540×256, 1920×1080×128, 3840×2160×64, and 1920×1080×32 respectively (the last deconvolution stride is 1, adjusting the size to the original 1920×1080). Finally, a 1×1 convolutional layer maps the number of channels from 32 to 3 (RGB three channels), generating an initial decoded image with a resolution of 1920×1080.

[0076] The detail enhancement algorithm employs a combined strategy of guided filtering and Laplacian enhancement. Guided filtering smooths the image while preserving edges, while Laplacian enhancement strengthens local texture details. The parameters for guided filtering are set as follows: window radius r = 5, regularization parameter ε = 0.01 (controlling the smoothing degree; a smaller ε results in a stronger smoothing effect). The initial decoded image is used as the guide map and filtered using the formula I_guide = GF(I_decode, I_decode, r, ε), where GF is the guided filtering function and I_decode is the initial decoded image. After filtering, the edge preservation rate is ≥95%, and the noise removal rate is ≥80%. In the example, the ground texture edges in the initial decoded image are blurred. After guided filtering, the edge gradient value increases from 25 to 38, significantly improving edge sharpness. Simultaneously, the overall noise grayscale value of the image decreases from 15 to 5, demonstrating good smoothing effect.

[0077] Laplacian enhancement extracts image details and textures by constructing a Laplacian operator. The enhancement formula is I_enhance = I_guide + k × L(I_guide), where L(I_guide) is the detail image after Laplacian operator processing, and k is the enhancement coefficient (set to 0.8 to avoid over-enhancement leading to noise amplification). The Laplacian operator uses a 3×3 convolution kernel with kernel values ​​of [0,-1,0;-1,4,-1;0,-1,0], which can effectively extract edge details of textures. In the example, the grayscale value of a certain region in the guided filter image I_guide is 145, while the corresponding region in the Laplacian detail image has a grayscale value of 8. After enhancement, I_enhance = 145 + 0.8 × 8 = 145 + 6.4 = 151.4, rounded to 152. The grayscale difference of texture details is increased from 8 to 13.4, resulting in clearer details.

[0078] The final image's illumination consistency is verified using a global illumination difference index. This calculates the difference in average grayscale values ​​between any two regions of the image; a difference ≤2% satisfies the global illumination consistency requirement. In the example, the average grayscale value of the highlight region in the final image is 232, the average grayscale value of the shadow region is 148, and the average grayscale value of the non-highlight and non-shadow regions is 150. Because the highlight region itself is brighter, the difference between the restored region and the non-highlight and non-shadow regions should be calculated as (150-148) / 148≈1.35%≤2%, satisfying global illumination consistency. The verification of local detail authenticity uses a texture similarity index. This compares the texture features (contrast and entropy of the grayscale co-occurrence matrix) of the final image and the original non-shadow region image. The final image has a contrast of 88 and an entropy of 1.92, while the original non-shadow region has a contrast of 85 and an entropy of 1.88. The difference ≤3.5% indicates that the local texture is realistic and reliable.

[0079] Finally, the image is normalized in grayscale and corrected in color gamut, limiting the grayscale range to [0, 255]. Gamma correction (Gamma value = 1.2) is used to adjust the image brightness to ensure a natural visual effect. The final output is an optical feature restoration image of the shadow area of ​​the land environment with global illumination consistency and local detail realism.

[0080] Another embodiment of the present invention provides a system for restoring optical features of shadowed areas in a terrestrial environment, see [link to documentation]. Figure 3 The system may include: The acquisition module 301 is used to simultaneously acquire visible light images, near-infrared spectral images and three-dimensional point cloud data of the target land area through a multimodal sensor, and generate multi-source fusion data with spatiotemporal alignment attributes; The construction module 302 is used to construct a physical model of light transmission in the scene based on the multi-source fusion data, and to separate the direct light attenuation component and the indirect light diffuse reflection component in the shadow area by solving the reflectivity equation under non-uniform light field. The reconstruction module 303 is used to compensate and reconstruct the direct illumination attenuation component using a generative adversarial network, and to perform spectral consistency correction based on the indirect illumination diffuse reflection component to generate an initial restored image after shadow removal. The output module 304 is used to adaptively weight and fuse the initial restored image with the original highlight region features using a cross-scale attention fusion mechanism, and output a final restored image with global illumination consistency and local detail authenticity.

[0081] This invention also provides an electronic device, including a memory and a processor, wherein the memory stores a computer program, and the processor is configured to run the computer program to perform the steps in any of the above method embodiments.

[0082] The above description, based on the embodiments shown in the figures, details the structure, features, and effects of the present invention. The above description is only a preferred embodiment of the present invention, but the present invention is not limited to the scope of implementation shown in the figures. Any changes made in accordance with the concept of the present invention, or equivalent embodiments modified to have equivalent changes, that do not exceed the spirit covered by the specification and figures, should be within the protection scope of the present invention.

Claims

1. A method for restoring optical features of shadowed areas in a terrestrial environment, characterized in that, The method includes: Visible light images, near-infrared spectral images, and 3D point cloud data of the target land area are acquired simultaneously by multimodal sensors to generate multi-source fusion data with spatiotemporal alignment attributes; Based on the multi-source fusion data, a physical model of light transmission in the scene is constructed. By solving the reflectivity equation under non-uniform illumination field, the direct light attenuation component and the indirect light diffuse reflection component in the shadow area are separated. Generative adversarial networks are used to compensate and reconstruct the direct illumination attenuation component, and spectral consistency correction is performed based on the indirect illumination diffuse reflection component to generate an initial restored image after shadow removal. A cross-scale attention fusion mechanism is adopted to adaptively weight and fuse the initial restored image with the original highlight region features, and output the final restored image with global illumination consistency and local detail authenticity.

2. The method according to claim 1, characterized in that, The process involves simultaneously acquiring visible light images, near-infrared spectral images, and 3D point cloud data of the target land area using multimodal sensors to generate multi-source fused data with spatiotemporal alignment attributes, including: A visible light camera, a near-infrared spectrometer, and a lidar sensor are integrated on the same observation platform. Hardware triggering circuits ensure that all sensors start acquiring data at the same time, generating a raw, synchronous multimodal data stream. The original synchronous multimodal data stream is timestamped and its spatial coordinate system is unified. The pre-calibrated sensor extrinsic parameter matrix is ​​used to transform the data of each mode to the same spatial coordinate system, generating multimodal data with unified spatiotemporal reference. Based on multimodal data with unified spatiotemporal reference, a feature point matching algorithm is used to accurately register 3D point clouds with visible light images and near-infrared spectral images, realizing a one-to-one correspondence between pixels and point clouds, and generating registered multi-source data; Multi-dimensional feature fusion is performed on the registered multi-source data to extract RGB channel features of visible light images, reflectance features of near-infrared spectral images, and geometric elevation features of three-dimensional point clouds, generating multi-source fused data with spatiotemporal alignment attributes.

3. The method according to claim 2, characterized in that, Based on the multi-source fused data, a physical model of light transmission in the scene is constructed. By solving the reflectivity equation under a non-uniform illumination field, the direct illumination attenuation component and the indirect illumination diffuse reflection component in the shadow region are separated, including: The three-dimensional geometric structure, surface normal vector distribution and material reflectivity characteristics of the scene are extracted from multi-source fusion data. Combined with solar azimuth angle, elevation angle and atmospheric transmittance parameters, a physical model of light transmission under non-uniform illumination field is constructed. Based on the physical model of light transmission, a reflectivity equation for the shadow region is established. This equation describes the attenuation process of direct light after it is blocked and the contribution of diffuse reflection formed by indirect light through multiple reflections. Using visible light images and near-infrared spectral images from multi-source fusion data as real observations, the least squares optimization algorithm is used to solve the reflectance equation to obtain the direct illumination attenuation coefficient and indirect illumination diffuse reflectance coefficient of the shadow region. Based on the attenuation coefficient and diffuse reflection coefficient obtained from the solution, the image data of the shadow area is decomposed into a direct illumination attenuation component and an indirect illumination diffuse reflection component, generating two independent illumination component images.

4. The method according to claim 3, characterized in that, The step of using a generative adversarial network to compensate and reconstruct the direct illumination attenuation component, and performing spectral consistency correction based on the indirect illumination diffuse reflection component to generate an initial restored image after shadow removal includes: A generative adversarial network model is constructed, in which the generator takes the direct illumination attenuation component as input and the discriminator takes the non-shaded area image as reference. The illumination compensation mapping relationship is learned through adversarial training to obtain the trained generative adversarial network. The direct illumination attenuation component is input into the trained generative adversarial network to generate the compensated and reconstructed direct illumination component, which restores the illumination intensity that should be in the shadow area. Using the diffuse reflection component of indirect illumination as the spectral reference, the spectral difference between the direct illumination component and the adjacent non-shaded area after compensation and reconstruction is calculated, and a spectral correction algorithm is used to eliminate spectral distortion in the compensation process. The spectrally corrected direct illumination component is linearly superimposed with the indirect illumination diffuse reflection component to generate an initial restored image after shadow removal, which has eliminated shadows but retains the original texture details.

5. The method according to claim 4, characterized in that, The method employs a cross-scale attention fusion mechanism to adaptively weight and fuse the initial restored image with the original highlight region features, outputting a final restored image with global illumination consistency and local detail fidelity, including: Multi-scale features of the highlight region are extracted from the original visible light image, and feature pyramids of the corresponding scale are extracted from the initial restored image to generate a highlight region feature set and an initial restored feature pyramid. Design a cross-scale attention fusion module that calculates the similarity between highlight region features and initial restored features at different scales and generates a scale-aware attention weight matrix. The initial restored feature pyramid is weighted and fused using an attention weight matrix, and the detailed features of the highlight region are adaptively fused into the initial restored image to generate a fused feature map. The fused feature map is decoded and reconstructed, the image resolution is restored through deconvolution, and the local texture is enhanced by a detail enhancement algorithm. The final output is a restored image with global illumination consistency and local detail realism.

6. A system for restoring optical features of shadowed areas in a terrestrial environment, characterized in that, The system includes: The acquisition module is used to simultaneously acquire visible light images, near-infrared spectral images, and three-dimensional point cloud data of the target land area through multimodal sensors, and generate multi-source fusion data with spatiotemporal alignment attributes; The construction module is used to construct a physical model of light transmission in the scene based on the multi-source fusion data, and to separate the direct light attenuation component and the indirect light diffuse reflection component in the shadow area by solving the reflectivity equation under non-uniform light field. The reconstruction module is used to compensate and reconstruct the direct illumination attenuation component using a generative adversarial network, and to perform spectral consistency correction based on the indirect illumination diffuse reflection component, thereby generating an initial restored image after shadow removal. The output module is used to adaptively weight and fuse the initial restored image with the original highlight region features using a cross-scale attention fusion mechanism, and output the final restored image with global illumination consistency and local detail authenticity.

7. The system according to claim 6, characterized in that, The acquisition module is specifically used for: A visible light camera, a near-infrared spectrometer, and a lidar sensor are integrated on the same observation platform. Hardware triggering circuits ensure that all sensors start acquiring data at the same time, generating a raw, synchronous multimodal data stream. The original synchronous multimodal data stream is timestamped and its spatial coordinate system is unified. The pre-calibrated sensor extrinsic parameter matrix is ​​used to transform the data of each mode to the same spatial coordinate system, generating multimodal data with unified spatiotemporal reference. Based on multimodal data with unified spatiotemporal reference, a feature point matching algorithm is used to accurately register 3D point clouds with visible light images and near-infrared spectral images, realizing a one-to-one correspondence between pixels and point clouds, and generating registered multi-source data; Multi-dimensional feature fusion is performed on the registered multi-source data to extract RGB channel features of visible light images, reflectance features of near-infrared spectral images, and geometric elevation features of three-dimensional point clouds, generating multi-source fused data with spatiotemporal alignment attributes.

8. The system according to claim 7, characterized in that, The building module is specifically used for: The three-dimensional geometric structure, surface normal vector distribution and material reflectivity characteristics of the scene are extracted from multi-source fusion data. Combined with solar azimuth angle, elevation angle and atmospheric transmittance parameters, a physical model of light transmission under non-uniform illumination field is constructed. Based on the physical model of light transmission, a reflectivity equation for the shadow region is established. This equation describes the attenuation process of direct light after it is blocked and the contribution of diffuse reflection formed by indirect light through multiple reflections. Using visible light images and near-infrared spectral images from multi-source fusion data as real observations, the least squares optimization algorithm is used to solve the reflectance equation to obtain the direct illumination attenuation coefficient and indirect illumination diffuse reflectance coefficient of the shadow region. Based on the attenuation coefficient and diffuse reflection coefficient obtained from the solution, the image data of the shadow area is decomposed into a direct illumination attenuation component and an indirect illumination diffuse reflection component, generating two independent illumination component images.

9. A storage medium, characterized in that, The storage medium stores a computer program, wherein the computer program is configured to execute the method of any one of claims 1-5 when it is run.

10. An electronic device comprising a memory and a processor, characterized in that, The memory stores a computer program, and the processor is configured to run the computer program to perform the method of any one of claims 1-5.

Images