Unmanned aerial vehicle remote sensing mapping image processing method based on image enhancement

Abstract

The application discloses an unmanned aerial vehicle remote sensing surveying and mapping image processing method based on image enhancement, and belongs to the technical field of unmanned aerial vehicle remote sensing surveying and mapping image processing. The method acquires an original remote sensing image sequence collected by a multispectral sensor carried by an unmanned aerial vehicle at multiple heading angles; generates an initial light direction field associated with the heading angle according to the radiation brightness value of each pixel point; extracts a terrain shadow strike feature by analyzing the vorticity component and the divergence component of the initial light direction field; performs radiation brightness reconstruction on the pixel points in the shadow coverage area according to the feature, and generates a shadow-removed remote sensing image; constructs a local texture enhancement weight of each pixel point based on the gradient modulus difference of adjacent pixel points in the shadow-removed image; and performs nonlinear stretching on the gray dynamic range of the shadow-removed image by using the weight, and outputs an enhanced remote sensing image.

Description

Technical Field

[0001] This invention relates to the field of UAV remote sensing mapping image processing technology, specifically to a UAV remote sensing mapping image processing method based on image enhancement. Background Technology

[0002] In UAV remote sensing mapping, multispectral sensors often acquire ground feature images at different heading angles along circular flight paths. Due to terrain undulations and changes in the sun's angle, shadows are prevalent in image sequences, leading to attenuated radiance and blurred textures in shadowed areas, thus reducing the accuracy of ground feature interpretation and mapping. Existing methods for shadow removal from remote sensing images largely rely on grayscale thresholding or classifier detection of shadow areas in single-frame images, followed by illumination recovery through histogram matching, brightness compensation, or Poisson editing. These approaches do not consider the systematic changes in illumination direction at multiple heading angles, making it difficult to accurately extract shadow directions related to terrain undulations. The recovered shadow areas often exhibit defects such as texture direction breaks, abrupt brightness changes, and unnatural transitions between non-shadowed areas. Some studies have attempted shadow reconstruction using digital surface models and illumination models, but these rely on high-precision exterior orientation elements and terrain data, resulting in complex processes and poor adaptability to low-overlap sequences. On the other hand, to improve the visual interpretation and feature extraction capabilities of remote sensing images, enhancement processing is usually required for the de-shadowed images. Existing enhancement techniques generally employ global linear stretching, histogram equalization, or gamma correction with fixed parameters, failing to adaptively adjust based on local texture differences. This can easily lead to the simultaneous amplification of noise in flat areas and insufficient contrast enhancement in high-texture areas, resulting in the obscuring of fine edges and subtle differences in ground features. In UAV multi-heading angle sequences, the radiance of the same ground feature pixel varies with the heading angle, forming a radiance vector. The first-order difference of this vector implies the relative relationship between the local illumination direction and the surface normal, which can be used to solve the horizontal projection of the illumination direction. However, directly obtaining the spatial extension features of shadows from the projection direction field is difficult, requiring a method that can stably describe the direction of terrain shadows. Simultaneously, how to use the extracted shadow direction features to guide the radiance interpolation of shadow areas, ensuring the natural continuity of the reconstructed texture with the surrounding terrain, is a problem that needs to be solved. Furthermore, for images after shadow removal, how to automatically adjust the stretching intensity based on the strength of local textures to enhance key details while avoiding the introduction of artifacts is also a problem in improving image quality. Summary of the Invention

[0003] This invention provides a UAV remote sensing mapping image processing method based on image enhancement to solve the problems in the prior art that do not make full use of multi-heading angle illumination information, resulting in texture distortion in the reconstruction of shadow areas, and the lack of local adaptive ability in the image enhancement process, leading to loss of detail.

[0004] To achieve the above objectives, the present invention provides the following technical solution: The present invention discloses an image processing method for UAV remote sensing mapping based on image enhancement, used to solve the problem of decreased accuracy of UAV remote sensing mapping caused by terrain shadows. The method includes: A sequence of raw remote sensing images is acquired from a multispectral sensor mounted on a UAV at multiple heading angles. An initial illumination direction field associated with the heading angle is generated based on the radiance values ​​of each pixel in the raw remote sensing image sequence. The curl and divergence components of the initial illumination direction field are analyzed to extract terrain shadow directional features from the raw remote sensing image sequence. Based on these terrain shadow directional features, radiance reconstruction is performed on pixels in the shadow-covered areas of the raw remote sensing image sequence to generate a de-shaded remote sensing image. Local texture enhancement weights are constructed for each pixel based on the gradient magnitude differences between adjacent pixels in the de-shaded remote sensing image. The grayscale dynamic range of the de-shaded remote sensing image is nonlinearly stretched using these local texture enhancement weights to output the enhanced remote sensing image.

[0005] Preferably, generating the initial illumination direction field includes: arranging the radiance values ​​of the same pixel at multiple heading angles into a radiance vector according to the heading angle; calculating the first-order difference of the radiance vector along the heading angle dimension to obtain the rate of change of radiance with the heading angle; and determining the projection angle of the illumination direction on the horizontal plane at each pixel based on the sign and magnitude of the rate of change, thus forming the initial illumination direction field. This method does not require a complex radiative transfer model, can efficiently and stably estimate the illumination direction, and provides a reliable foundation for subsequent shadow processing.

[0006] Furthermore, when the rate of change is positive, the projection angle is assigned the heading angle plus 90 degrees; when the rate of change is negative, it is assigned the heading angle minus 90 degrees. The amplitude is normalized and multiplied by a preset angle perturbation threshold to obtain the angle perturbation offset, which is then used to correct the projection angle. The corrected projection angles of each pixel are combined to form the initial illumination direction field. This correction enhances the robustness of the direction field to sensor noise and changes in ground reflectivity.

[0007] The process of extracting terrain shadow direction features is as follows: The lateral and longitudinal partial derivatives of the initial illumination direction field are calculated in the image plane coordinate system. The curl component is calculated based on the difference between the lateral and longitudinal partial derivatives, and the divergence component is calculated by summing the two. The curl component and the divergence component are multiplied to obtain a product matrix. Local maxima detection is performed on the product matrix, and the spatial distribution lines corresponding to the local maxima points represent the terrain shadow direction features. Using curl and divergence information, the directionality of shadows caused by terrain undulations can be accurately characterized. Compared with the simple brightness thresholding method, this significantly reduces false detections of low-reflectivity features such as vegetation and soil.

[0008] The de-shaded remote sensing image is generated by reconstructing radiance based on the terrain shadow directional features. Preferably, along the normal direction of the terrain shadow directional features, a non-shadow reference pixel closest to the shadow-covered pixel in the original remote sensing image sequence is searched. Based on the radiance value of the non-shadow reference pixel and its spatial distance to the shadow-covered pixel, an inverse distance interpolation is used to calculate the reconstructed radiance value, which replaces the original radiance value, thus obtaining the de-shaded remote sensing image. This reconstruction strategy follows the spatial gradation law of shadows, perfectly preserving terrain texture details while eliminating shadows.

[0009] As a preferred search method, two search rays are emitted from the center of the shadow-covered pixel, one in the forward direction and one in the reverse direction of the normal. The search ray expands outward pixel by pixel along each ray until a pixel with a radiance value greater than a preset shadow determination threshold is encountered; this pixel is then selected as a candidate reference pixel. The Euclidean distances between the candidate reference pixels found by the two search rays and the pixels in the shadow-covered area are compared, and the candidate reference pixel with the smaller distance is selected as the non-shadow reference pixel. This ensures that the selected reference pixels are most similar in lighting conditions, avoiding blurring caused by interpolation across texture boundaries.

[0010] When constructing local texture enhancement weights, preferably, the gradient components of each pixel in the deshaded remote sensing image in the horizontal and vertical directions are calculated, and the gradient magnitude is obtained; with each target pixel as the center, the average gradient magnitude of the neighboring pixels (excluding the target pixel) within a neighborhood window of a preset radius is calculated; the ratio of the gradient magnitude of the target pixel to the average gradient magnitude is used as the local texture enhancement weight for that pixel. This weight can significantly highlight the edges and textures of ground features, while suppressing the gain in smooth areas and preventing background noise from being excessively amplified.

[0011] Furthermore, the local texture enhancement weights are cropped: the weight values ​​of each pixel are sorted from smallest to largest, and the top 5% and bottom 5% of extreme weights are truncated and replaced with the minimum and maximum values ​​of the remaining weights, respectively, to obtain the cropped local texture enhancement weights. Non-linear stretching is then performed using these cropped weights. This process reduces the influence of abnormal extreme values, resulting in a more uniform and natural overall enhancement effect.

[0012] When performing nonlinear stretching using local texture enhancement weights, the preferred method is to divide the original grayscale value of each pixel in the shaded remote sensing image by the maximum grayscale value of the entire image to obtain a normalized grayscale value; then, using the local texture enhancement weights as the exponent, perform a power transformation on the normalized grayscale value to obtain an enhanced grayscale value; finally, multiply the enhanced grayscale value by the maximum grayscale value of the entire image and round it to obtain the output grayscale value, which constitutes the enhanced remote sensing image. Adaptive power transformation effectively stretches the contrast of dark areas such as shadows while maintaining the detail of bright areas, significantly improving image interpretability and mapping accuracy.

[0013] Regarding the acquisition of the original remote sensing image sequence, preferably, the UAV is controlled to fly along a preset circular track, triggering the multispectral sensor to simultaneously acquire visible light and near-infrared images at equally spaced angular positions; the visible light and near-infrared images acquired at the same angular position are then registered to form a multi-band original image; these images are arranged in ascending order of heading angle to form the original remote sensing image sequence. The circular track and multispectral registration provide complementary information from multiple angles and bands for the illumination direction field calculation, which is beneficial to improving the accuracy of shadow detection and removal.

[0014] The technical effects and advantages provided by the present invention in the above technical solution are as follows: A luminance vector is constructed based on the radiance values ​​of the same pixel at multiple heading angles. The first-order difference of this vector along the heading angle dimension is calculated to obtain the rate of change of radiance with the heading angle. The projection angle of the illumination direction onto the horizontal plane is assigned to the heading angle plus or minus 90 degrees based on the sign of the rate of change. The amplitude of the rate of change is then normalized and multiplied by a preset threshold to obtain the angle perturbation offset. After correction, an initial illumination direction field is formed. The lateral and longitudinal partial derivatives of this direction field are calculated. The curl component is obtained from the difference between the lateral and longitudinal partial derivatives, and the divergence component is obtained by summing the two. The curl component and the divergence component are multiplied to obtain a product matrix. Local maxima are detected in the product matrix, and the spatial distribution of these maxima points constitutes the terrain shadow direction feature. This method utilizes the luminance gradient information brought about by multi-angle illumination changes and directly extracts the spatial direction of terrain shadows through curl and divergence decomposition in field theory. It does not rely on external terrain elevation data, making the description of shadow direction closely related to the terrain undulation, exhibiting strong noise resistance and directional consistency. During the reconstruction phase, search rays are emitted from the shadowed pixels in both the forward and reverse directions along the normal direction of the shadow's direction feature. This process expands pixel by pixel until a non-shadowed reference pixel with a radiance value higher than a preset threshold is encountered. The pixel with the closest Euclidean distance is selected as the reference point. The reconstructed radiance value is calculated using inverse distance weighting and replaces the original value, generating a de-shadowed remote sensing image. Thanks to accurate shadow direction constraints, the search path is always perpendicular to the shadow stretching direction. The probability of non-shadowed reference points and the pixel to be restored being within the same terrain texture unit is higher, ensuring the interpolation results maintain the continuity of the texture direction. This avoids texture misalignment and brightness abrupt changes caused by improper reference point selection in traditional methods. The image after shadow removal has a smooth brightness transition and realistic ground feature morphology. The gradient components in the horizontal and vertical directions of the de-shadowed remote sensing image are calculated, and the gradient magnitude is obtained. Centered on each target pixel, the average gradient magnitude of neighboring pixels within a preset neighborhood window (excluding the target pixel) is calculated. The ratio of the target pixel's gradient magnitude to the average gradient magnitude is used as the local texture enhancement weight for that point. A ratio greater than 1 indicates that the texture of the pixel's region is more significant than its neighbors, while a ratio less than 1 indicates relative flatness. After sorting all weights, the top 5% and bottom 5% extreme values ​​are removed to avoid excessive influence of isolated points or abrupt noise on the stretching intensity. Using the cropped local texture enhancement weights as a power exponent, a power transformation is performed on the normalized grayscale values ​​of the deshaded image. The result is then mapped back to the original grayscale range to output the enhanced remote sensing image. Regions with rich texture have larger weights; as the power exponent increases, low grayscale areas are compressed and high grayscale areas are enhanced, significantly improving detail contrast. Flat regions have smaller weights, resulting in a near-linear grayscale transformation and suppressing noise amplification. This adaptive nonlinear stretching process is entirely driven by the image's own gradient structure, requiring no preset global parameters. This allows the enhanced image to possess both sharp edges in high-texture regions and soft layers in flat regions, improving the presentation of internal differences in ground features and facilitating subsequent ground feature classification and boundary delineation. Attached Figure Description

[0015] To more clearly illustrate the technical solutions in the embodiments of this application or the prior art, the drawings used in the embodiments will be briefly introduced below. Obviously, the drawings described below are only some embodiments recorded in this invention. For those skilled in the art, other drawings can be obtained based on these drawings.

[0016] Figure 1 This is a flowchart of an image processing method for UAV remote sensing mapping based on image enhancement; Figure 2 This is a flowchart of the UAV circular flight path multispectral image acquisition and registration process; Figure 3 This is a schematic diagram of the UAV's circular flight path and data collection trigger points; Figure 4 These are the curves showing the numerical changes of the curl component, divergence component, and product matrix. Detailed Implementation

[0017] To make the objectives, technical solutions, and advantages of the embodiments of the present invention clearer, the technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0018] See Figure 1 This invention provides an image processing method for UAV remote sensing mapping based on image enhancement. The method involves acquiring a sequence of original remote sensing images collected by a multispectral sensor mounted on a UAV at multiple heading angles; generating an initial illumination direction field associated with the heading angle based on the radiance values ​​of each pixel in the original remote sensing image sequence; extracting terrain shadow direction features from the original remote sensing image sequence by analyzing the curl and divergence components of the initial illumination direction field; reconstructing the radiance of pixels in the shadow-covered area of ​​the original remote sensing image sequence based on the terrain shadow direction features to generate a shaded remote sensing image; constructing local texture enhancement weights for each pixel based on the gradient magnitude differences between adjacent pixels in the shaded remote sensing image; and using these local texture enhancement weights to nonlinearly stretch the grayscale dynamic range of the shaded remote sensing image to output the enhanced remote sensing image. Example

[0019] In specific implementation, please refer to Figure 2The system controls the UAV to fly along a pre-set circular trajectory within the target survey area. This trajectory is generated by the ground control station based on the geographical location and extent of the target survey area. The radius of the circular trajectory is set according to the field of view of the multispectral sensor and the required ground sampling distance. The UAV flight control system loads the waypoint coordinates of the circular trajectory and adjusts the UAV's flight attitude and position in real time via satellite navigation positioning, ensuring the UAV flies at a constant speed along the circular trajectory. During flight, the UAV's heading angle changes continuously with the circular trajectory.

[0020] At equally spaced angular positions along the circular flight path, the multispectral sensor is triggered to simultaneously acquire images in the visible and near-infrared bands. The equally spaced angles are determined based on the total change in the heading angle and the preset number of image acquisitions. The angle interval is calculated by dividing 360 degrees by the preset number of acquisitions. The flight control system monitors the UAV's current heading angle in real time. Whenever the change in the current heading angle reaches an integer multiple of the angle interval value, a trigger signal is generated and sent to the multispectral sensor. Upon receiving the trigger signal, the multispectral sensor simultaneously activates the visible and near-infrared channels at the same exposure time to acquire visible and near-infrared images of the same ground feature scene. Each triggered acquisition corresponds to a specific acquisition position, and the heading angle value corresponding to that acquisition position is recorded.

[0021] In practice, visible light and near-infrared images acquired at the same angular location are registered together. The registration process involves using the visible light image as a reference image and aligning the near-infrared image to the pixel grid of the visible light image through geometric transformation. The steps for obtaining the geometric transformation parameters include: extracting several pairs of corresponding feature points from both the reference image and the near-infrared image to be registered; matching these pairs using a scale-invariant feature transformation algorithm; solving for the affine transformation matrix using the matched pairs of feature points, which includes translation, rotation, and scaling parameters; and resampling the near-infrared image according to the affine transformation matrix, ensuring that each pixel in the near-infrared image represents the same ground location as its corresponding pixel in the visible light image. After registration, the visible light image generated in the same trigger and the registered near-infrared image are merged at the band level to form a multi-band original image for that angular location. Each pixel in the multi-band original image contains at least a visible light radiance value and a near-infrared radiance value.

[0022] In practice, all multi-band raw images at all angular positions are arranged in ascending order of heading angle to form a raw remote sensing image sequence. The flight control system associates the heading angle value recorded each time acquisition is triggered with the corresponding multi-band raw image. All multi-band raw images are then sorted in ascending order of their associated heading angle values ​​to obtain the raw remote sensing image sequence. Each frame of the multi-band raw image in the raw remote sensing image sequence corresponds to a known heading angle.

[0023] See Figure 3 The figure shows the distribution of the UAV's pre-set circular flight path and actual flight trajectory points within the target survey area. The horizontal axis represents the lateral X-coordinate (unit: meters), and the vertical axis represents the longitudinal Y-coordinate (unit: meters), both ranging from approximately -1100 meters to 1100 meters. The dashed line in the figure represents the pre-set circular flight path, which is approximately a perfect circle and delineates the ideal flight path of the UAV. The scattered points, represented by small blue dots, represent the track points collected during the actual flight of the UAV. The actual track points closely surround the pre-set circular flight path, indicating that the UAV can fly at a constant speed along the predetermined path relatively well.

[0024] The actual acquisition trigger points are marked with brown-red squares in the figure. These trigger points are evenly distributed along the circular track, and their number and location correspond to the equally spaced acquisition angles of the preset circular track, verifying the real-time generation mechanism of the acquisition trigger signal. The uniform intervals of the trigger points indicate that the flight control system accurately triggers multispectral sensor acquisition based on changes in heading angle, ensuring the continuity and uniformity of heading angle sampling for each frame in the multi-band original image sequence. Example

[0025] In practice, the radiance values ​​of the same pixel at multiple heading angles are extracted from the original remote sensing image sequence. For the image plane coordinate system with coordinates... For each pixel in the original remote sensing image sequence, the radiance value of that pixel in the visible light band is obtained by traversing each frame of the original multi-band image. The obtained radiance values ​​are then arranged in ascending order according to the heading angle corresponding to the multi-band image, forming a radiance vector. The radiance vector is represented as follows: ,in Indicates coordinates as The pixel in the first a heading angle The corresponding radiance value, This indicates the total number of multi-band original image frames contained in the original remote sensing image sequence. The value is to Integers.

[0026] In practice, the first-order difference of the brightness vector along the heading angle dimension is calculated to obtain the rate of change of radiance with the heading angle. For coordinates... For each pixel, the difference in radiance value between adjacent heading angles is calculated using the following formula:

[0027] in, Indicates coordinates as The pixel in the first a heading angle The rate of change of radiance at that location, Indicates the first One heading angle, Indicates coordinates as The pixel in the first a heading angle The corresponding radiance value. For from Get the value to Calculated sequentially The radiance change rate constitutes a radiance change rate sequence.

[0028] In practice, the projection angle of the illumination direction on the horizontal plane at each pixel is determined based on the sign and magnitude of the rate of change. For coordinates... For each pixel in the radiance change rate sequence, traverse each radiance change rate. When the rate of change of radiance When the sign is positive, the initial projection angle is assigned the heading angle. Add ninety degrees; when the rate of change of radiance When the sign of the value is negative, the initial projection angle is assigned the heading angle. Reduce by 90 degrees.

[0029] In practice, the amplitude of the rate of change of radiance is normalized. All rates of change of radiance in the radiance change rate sequence are then obtained. Find the absolute value of the given information and determine its maximum value. and minimum value For each rate of change of radiance The amplitude is normalized, and the normalized amplitude is... The calculation method is as follows: The normalized amplitude ranges from zero to one.

[0030] In practice, a preset angle disturbance threshold is obtained. This threshold is set based on equally spaced angle values ​​of the heading angle, and is set to half of each equally spaced angle value. The equally spaced angle value is calculated by dividing 360 degrees by the total number of frames in the original remote sensing image sequence. We obtain the normalized amplitude. Multiply by the angle perturbation threshold to obtain the angle perturbation offset.

[0031] In practice, the initial projection angle is added to the angle disturbance offset to obtain the corrected projection angle. For each heading angle... The corrected projection angle is obtained by adding the corresponding initial projection angle to the corresponding angle perturbation offset. All heading angles Corresponding corrected projection angle The coordinates are obtained by averaging. The final projection angle of the pixel The final projection angle is calculated for each pixel in the original remote sensing image sequence, and the final projection angles of each pixel are combined according to their positions in the image plane coordinate system to form the initial illumination direction field. Example

[0032] In practice, the lateral and longitudinal partial derivatives of the initial illumination orientation field are calculated in the image plane coordinate system. The initial illumination orientation field is stored as a two-dimensional matrix, where the row coordinates correspond to the longitudinal coordinates of the image plane coordinate system, the column coordinates correspond to the lateral coordinates, and the value of each element in the matrix represents the final projection angle of the corresponding pixel. The lateral partial derivative is calculated as follows: for the initial illumination orientation field with coordinates... For a given pixel, subtract the final projection angle of its right-hand neighbor from the final projection angle of its left-hand neighbor, and divide the difference by twice the pixel spacing to obtain its coordinates. The horizontal partial derivative of the pixel. The vertical partial derivative is obtained as follows: for the initial illumination direction field with coordinates... For a given pixel, subtract the final projection angle of its adjacent pixel below from the final projection angle of its adjacent pixel above. Divide the difference by twice the pixel spacing to obtain the coordinates. The vertical partial derivatives of the pixels are calculated. The pixel spacing is the actual geographic distance between the centers of adjacent pixels in the original multi-band image. For pixels at the initial illumination direction field boundary, the lateral and vertical partial derivatives are calculated using a one-sided difference method.

[0033] In practical implementation, the curl component is calculated based on the difference between the transverse and longitudinal partial derivatives. For coordinates... Pixels, curl components The calculation method is shown in the following formula:

[0034] in, Indicates coordinates as The curl component of the pixel, The coordinates in the initial illumination direction field are: The horizontal component of the final projection angle of the pixel. Indicates the horizontal component along the vertical coordinate. Partial derivatives in direction, The coordinates in the initial illumination direction field are: The vertical component of the final projection angle of the pixel. Indicates the vertical component along the horizontal coordinate. Partial derivative in direction. Lateral component. The value is the cosine of the final projected angle, and the longitudinal component is... The value of is the sine of the final projection angle.

[0035] In practical implementation, the divergence component is calculated by summing the horizontal and vertical partial derivatives. For coordinates... pixel, divergence component The calculation method is as follows: the horizontal component For the horizontal coordinate Find the partial derivatives to obtain the first partial derivative; then convert the vertical components... For the vertical axis Find the partial derivative to obtain the second partial derivative; add the first and second partial derivatives to obtain the divergence components. The first partial derivative is obtained by taking the coordinates as... The second partial derivative is calculated by subtracting the horizontal component of the right-hand neighboring pixel from the horizontal component of the left-hand neighboring pixel, and then dividing the difference by twice the pixel spacing. The second partial derivative is obtained by taking the coordinates as... The vertical component of the pixel below the pixel is subtracted from the vertical component of the pixel above the pixel, and the difference is divided by twice the pixel spacing.

[0036] In practice, the curl component is multiplied by the divergence component to obtain the product matrix. For each coordinate... Calculate the product value of the pixels. , product value The value is assigned to the element at the corresponding position in the product matrix. The product value is calculated for each pixel in the initial illumination direction field to construct the complete product matrix.

[0037] In practice, local maxima detection is performed on the product matrix. The process of local maxima detection is as follows: for a product matrix with coordinates... The elements, in the coordinates of Within a preset-sized neighborhood window centered on the element, the product value of that element is compared with the product values ​​of all elements within the neighborhood window. The preset-sized neighborhood window has a size of three rows and three columns, meaning it contains elements at coordinates of... The elements and their adjacent elements in the top, bottom, left, right, and four diagonal directions. If the coordinates are... If the product of the elements is greater than the product of all adjacent elements within the neighborhood window, then the coordinates are... The elements are marked as local maxima. The local maxima are detected by traversing all elements of the product matrix, resulting in a set of local maxima.

[0038] In practice, the spatial distribution lines corresponding to local maxima are used as the terrain shadow direction features. Each local maximum point in the set of local maximum points corresponds to a spatial location in the image plane coordinate system. All local maximum points are connected according to their spatial adjacency. The connection method is as follows: for each local maximum point, search whether there are other local maximum points within its eight-neighborhood. If so, connect the local maximum point with the local maximum points in its neighborhood using line segments. All line segments combine to form multiple continuous spatial distribution lines. The direction and extension range of these spatial distribution lines in the image plane coordinate system constitute the terrain shadow direction features.

[0039] See Figure 4 In the graph, the horizontal axis represents the pixel index, ranging from approximately 0 to 600, and the vertical axis represents the curl component. , divergence component and their product The dimensionless value. The green solid line curve represents the curl component. The purple dashed curve represents the divergence component. The red dotted line curve represents the product matrix of the two. The curl component curve exhibits significant periodic fluctuations within the pixel index range, with values ​​varying between approximately ±0.7, reflecting the rotational characteristics of the initial illumination orientation field on the image plane. The divergence component curve also shows periodic variations, with slightly smaller amplitudes, ranging from approximately ±0.55, characterizing the divergence or convergence of the illumination orientation field. The product matrix curve shows significantly lower fluctuation amplitudes than the curl and divergence components, with values ​​concentrated within the ±0.1 range, and exhibits the characteristic of crossing the zero axis multiple times, indicating local variations in the curl-divergence product. Example

[0040] In the specific implementation, the direction of each spatially distributed line in the terrain shadow direction feature is obtained. For each line segment on the spatially distributed line constituting the terrain shadow direction feature, the direction vector of the line segment is calculated. The direction vector is calculated by subtracting the coordinates of the starting point of the line segment from the coordinates of the ending point. The normal direction is calculated based on the direction vector. The normal direction is perpendicular to the direction vector. The normal direction is calculated by rotating the direction vector counterclockwise by 90 degrees in the image plane coordinate system to obtain the positive normal direction vector, and rotating the direction vector clockwise by 90 degrees in the image plane coordinate system to obtain the negative normal direction vector. The positive and negative normal direction vectors are then normalized to obtain the positive unit normal vector and the negative unit normal vector, respectively.

[0041] In practice, the pixels in the shadowed areas of the original remote sensing image sequence are identified. The method for determining shadowed pixels is as follows: For each frame of the multi-band original image in the original remote sensing image sequence, the average visible light radiance value of all pixels in the entire image is calculated. This average radiance value is then multiplied by a preset shadow determination coefficient to obtain a preset shadow determination threshold. The preset shadow determination coefficient is set to 0.3. The reason for this value is that in terrain shadow areas, the direct solar radiation received by the surface is completely blocked. The surface radiance in shadow areas mainly comes from atmospheric scattering and reflected radiation from surrounding objects. Its radiance value typically drops to 20% to 40% of the radiance value in directly sunny areas. Using 0.3 as the shadow determination coefficient corresponds to the typical radiance boundary ratio between shadowed and non-shadowed areas. Each pixel in the multi-band original image is iterated over; if the visible light radiance value of a pixel is less than the preset shadow determination threshold, then that pixel is marked as a shadowed pixel.

[0042] In practice, for each pixel in the shadow-covered area, a non-shadow reference pixel is searched along the normal direction of the terrain shadow feature. Centered on the pixel in the shadow-covered area, a forward search ray is emitted along the positive unit normal vector direction, and a reverse search ray is emitted along the negative unit normal vector direction. The forward search ray is generated as follows: starting from the coordinates of the pixel in the shadow-covered area, it moves pixel by pixel along the positive unit normal vector direction. The coordinates of the pixel reached in each step are obtained by adding the coordinates of the previous step to the positive unit normal vector and rounding down. The reverse search ray is generated as follows: starting from the coordinates of the pixel in the shadow-covered area, it moves pixel by pixel along the negative unit normal vector direction. The coordinates of the pixel reached in each step are obtained by adding the coordinates of the previous step to the negative unit normal vector and rounding down.

[0043] In practice, the search for candidate reference pixels expands outward pixel by pixel along each search ray. After each step to a new pixel, the forward search ray acquires the visible light radiance value of that pixel and compares it with a preset shadow determination threshold. If the visible light radiance value is greater than the threshold, the forward search ray stops expanding and the pixel is designated as a forward candidate reference pixel. If the visible light radiance value is less than or equal to the threshold, the ray continues to move along the forward unit normal vector to the next pixel, repeating the radiance value comparison process. The reverse search ray employs the same search strategy: after each step to a new pixel, the visible light radiance value is acquired and compared with a preset shadow determination threshold. If the value is greater than the threshold, the ray stops expanding and the pixel is designated as a reverse candidate reference pixel; otherwise, the ray continues to move forward.

[0044] In practice, the Euclidean distances between the forward and reverse candidate reference pixels and the shadow-covered pixels are compared. The Euclidean distance between the forward candidate reference pixel and the shadow-covered pixel is calculated as follows: the square of the difference between the row coordinates of the forward candidate reference pixel and the row coordinates of the shadow-covered pixel is added to the square of the difference between the column coordinates of the forward candidate reference pixel and the column coordinates of the shadow-covered pixel, and the square root of the sum is taken to obtain the forward Euclidean distance. Similarly, the Euclidean distance between the reverse candidate reference pixel and the shadow-covered pixel is calculated as follows: the square of the difference between the row coordinates of the reverse candidate reference pixel and the row coordinates of the shadow-covered pixel is added to the square of the difference between the column coordinates of the reverse candidate reference pixel and the column coordinates of the shadow-covered pixel, and the square root of the sum is taken to obtain the reverse Euclidean distance. The forward and reverse Euclidean distances are compared, and the candidate reference pixel with the smaller Euclidean distance is selected as the non-shadow reference pixel corresponding to the shadow-covered pixel.

[0045] In practice, based on the radiance value of the non-shaded reference pixel and the spatial distance between the shadowed pixels and the non-shaded reference pixels, the reconstructed radiance value of the shadowed pixels is calculated using inverse distance interpolation. The calculation method is shown in the following formula:

[0046] in, Indicates coordinates as The reconstructed radiance value of the pixels in the shadow-covered area. This represents the total number of non-shadow reference pixels involved in the interpolation calculation. The value of is 1, meaning that interpolation is performed using only the nearest neighbor single non-shaded reference pixel. Indicates the first The non-shadow reference pixels are at coordinates... The visible light radiance value at that location. Indicates the pixel of the shadow coverage area and the first pixel. The spatial distance between non-shadow reference pixels is calculated using Euclidean distance.

[0047] In practice, the reconstructed radiance value replaces the original radiance value of the corresponding pixel in the original remote sensing image sequence. For each pixel in the shadowed area, the calculated reconstructed radiance value is... Assigning the coordinates to the original multi-band image The visible light radiance values ​​of the pixels in the original multi-band image are replaced one by one, overwriting the original visible light radiance values. After replacing the radiance values ​​of all pixels in the shadowed areas of the original multi-band image, a shadow-free remote sensing image is obtained. Pixels in the non-shadowed areas of the shadow-free remote sensing image retain their original visible light radiance values. Example

[0048] In practice, the gradient components of each pixel in the deshaded remote sensing image are calculated in the horizontal and vertical directions. The horizontal gradient component is calculated as follows: for a pixel in the deshaded remote sensing image with coordinates... For a given pixel, subtract the grayscale value of its right-hand neighbor from the grayscale value of its left-hand neighbor, and divide the difference by two to obtain its coordinates. Horizontal gradient components of pixels The vertical gradient component is calculated as follows: For a shaded remote sensing image with coordinates... For a given pixel, subtract the grayscale value of the adjacent pixel below it from the grayscale value of the adjacent pixel above it, and divide the difference by two to obtain the coordinates. Vertical gradient components of pixels For pixels at the boundaries of the shadow-removed remote sensing image, the horizontal and vertical gradient components are calculated using a one-sided difference method. Specifically, for pixels at the left boundary of the image, the gray value of the current pixel is taken by subtracting the gray value of the adjacent pixel on the right as the horizontal gradient component. For pixels at the right boundary of the image, the gray value of the current pixel is taken by subtracting the gray value of the adjacent pixel on the left as the horizontal gradient component. The boundary processing method for the vertical gradient component is similar.

[0049] In practice, the gradient magnitude of each pixel is calculated based on the horizontal and vertical gradient components. For a shadow-removed remote sensing image with coordinates... Pixels, gradient magnitude The calculation method is as follows: calculate the gradient component in the horizontal direction. The square of the gradient component in the vertical direction The sum of squares, taking the square root of the sum, yields the coordinates. The gradient magnitude of each pixel is calculated. This is done by iterating through all pixels in the deshaded remote sensing image and calculating the gradient magnitude of each pixel to obtain the gradient magnitude image.

[0050] In practice, for each target pixel, the average gradient magnitude of all neighboring pixels (excluding the target pixel) within a neighborhood window of a preset radius is calculated. The preset radius is set to three pixels. This is because a three-pixel radius neighborhood window can cover a sufficient number of neighboring pixels around the target pixel to obtain a statistically stable average gradient magnitude, while avoiding an excessively large neighborhood window that would cause a smooth transition of statistical values ​​across different texture regions. The neighborhood window is a seven-row, seven-column square area, extending three pixels upwards, downwards, leftwards, and rightwards from the target pixel. For each pixel within the neighborhood window, if its coordinates differ from the target pixel's coordinates, it is included in the neighborhood pixel set. The gradient magnitudes of all neighboring pixels in the neighborhood pixel set are calculated, and the sum of these gradient magnitudes is divided by the number of pixels in the neighborhood pixel set to obtain the average gradient magnitude. .

[0051] In practice, the ratio of the gradient magnitude of the target pixel to its average gradient magnitude is used as the local texture enhancement weight for that pixel. The calculation method is shown in the following formula:

[0052] in, Indicates coordinates as Local texture enhancement weights for target pixels. Indicates coordinates as The gradient magnitude of the target pixel. Represented by coordinates The average gradient magnitude of all neighboring pixels (excluding the target pixel) within a preset radius neighborhood window centered on the target pixel is calculated. Local texture enhancement weights are then calculated for each pixel in the deshaded remote sensing image to obtain a local texture enhancement weight image.

[0053] In practice, the local texture enhancement weights of each pixel are sorted from smallest to largest. The local texture enhancement weight values ​​of all pixels in the local texture enhancement weight image are obtained, forming a weight value sequence. The length of the weight value sequence is equal to the total number of pixels in the shadow-removed remote sensing image. A sorting algorithm is used to arrange all values ​​in the weight value sequence in ascending order, resulting in a sorted weight sequence.

[0054] In practice, the extreme weights located in the first 5% and the last 5% of the sorted weight sequence are extracted. The extraction method for the first 5% extreme weights is as follows: calculate the extraction quantity, which is equal to 5% of the total length of the sorted weight sequence and rounded down. Starting from the beginning of the sorted weight sequence, extract a certain number of weight values ​​as the first 5% extreme weights. The extraction method for the last 5% extreme weights is as follows: starting from the end of the sorted weight sequence, extract a certain number of weight values ​​backwards as the last 5% extreme weights.

[0055] In practice, the extreme weights are replaced with the minimum and maximum values ​​from the remaining weights. The remaining weights refer to the middle weight values ​​retained after removing the top 5% and bottom 5% extreme weights from the sorted weight sequence. The minimum value is found from the remaining weights, and every value in the top 5% extreme weights is replaced with the minimum value from the remaining weights. Similarly, the maximum value is found from the remaining weights, and every value in the bottom 5% extreme weights is replaced with the maximum value from the remaining weights. After the replacements are complete, the sorted weight sequence is restored to its corresponding position in the original image's planar coordinate system, resulting in the cropped local texture enhancement weighted image.

[0056] In practice, the cropped local texture enhancement weights are used to replace the original local texture enhancement weights for nonlinear stretching. The normalized grayscale value is obtained by dividing the original grayscale value of each pixel in the deshaded remote sensing image by the maximum grayscale value of the entire image. The maximum grayscale value of the entire image is obtained by iterating through all pixels in the deshaded remote sensing image and finding the maximum grayscale value. For a pixel in the deshaded remote sensing image with coordinates... The pixels, normalized gray values The calculation method is as follows: The coordinates are... The original grayscale value of each pixel is divided by the maximum grayscale value of the entire image.

[0057] In practice, the local texture enhancement weights after cropping are used as the exponent to perform a power-law transformation on the normalized grayscale values ​​to obtain the enhanced grayscale values. For the deshading remote sensing image with coordinates... Pixels, enhance grayscale values The calculation method is as follows: normalize the gray value The value is exponentially calculated, with the exponent being the coordinates. The local texture enhancement weight corresponding to the pixel after cropping The mathematical expression for the power transformation is: of Power of 1.

[0058] In practice, the enhanced grayscale value is multiplied by the maximum grayscale value of the entire image and then rounded to obtain the output grayscale value. For a shadow-removed remote sensing image with coordinates... The pixels will enhance the grayscale value. Multiply by the maximum grayscale value of the entire image, round the product to the nearest integer, and obtain the coordinates. The output grayscale value of each pixel is calculated. The enhanced remote sensing image is then constructed by combining the output grayscale values ​​of all pixels according to their corresponding positions in the image's planar coordinate system.

[0059] The above description is merely a specific embodiment of this application, but the scope of protection of this application is not limited thereto. Any changes or substitutions that can be easily conceived by those skilled in the art within the scope of the technology disclosed in this application should be included within the scope of protection of this application.

Claims

1. A method for processing UAV remote sensing and mapping images based on image enhancement, characterized in that, include: Acquire raw remote sensing image sequences from the multispectral sensor onboard the UAV at multiple heading angles; Based on the radiance values ​​of each pixel in the original remote sensing image sequence, an initial illumination direction field associated with the heading angle is generated; By analyzing the curl and divergence components of the initial illumination direction field, the terrain shadow orientation features in the original remote sensing image sequence are extracted. Based on the terrain shadow direction features, the radiance reconstruction is performed on the pixels in the shadow-covered area of ​​the original remote sensing image sequence to generate a shadow-free remote sensing image. Based on the difference in gradient magnitude between adjacent pixels in the deshaded remote sensing image, local texture enhancement weights are constructed for each pixel. The grayscale dynamic range of the deshaded remote sensing image is nonlinearly stretched using the local texture enhancement weights to output the enhanced remote sensing image.

2. The UAV remote sensing mapping image processing method based on image enhancement according to claim 1, characterized in that, The step of generating an initial illumination direction field associated with the heading angle based on the radiance values ​​of each pixel in the original remote sensing image sequence includes: Arrange the radiance values ​​of the same pixel in the original remote sensing image sequence at multiple heading angles into a radiance vector according to the heading angle order; Calculate the first-order difference of the brightness vector in the heading angle dimension to obtain the rate of change of radiance with heading angle; Based on the sign and magnitude of the rate of change, the projection angle of the illumination direction on the horizontal plane at each pixel is determined, thus forming the initial illumination direction field.

3. The UAV remote sensing and mapping image processing method based on image enhancement according to claim 1, characterized in that, The step of extracting terrain shadow orientation features from the original remote sensing image sequence by analyzing the curl and divergence components of the initial illumination direction field includes: The lateral and longitudinal partial derivatives of the initial illumination direction field are obtained in the image plane coordinate system, and the curl component is calculated based on the difference between the lateral and longitudinal partial derivatives. The divergence component is calculated by summing the lateral and longitudinal partial derivatives. The curl component and the divergence component are multiplied to obtain a product matrix, and local maxima detection is performed on the product matrix. The spatial distribution lines corresponding to the local maxima points are used as the terrain shadow direction features.

4. The UAV remote sensing and mapping image processing method based on image enhancement according to claim 1, characterized in that, The step of reconstructing the radiometric brightness of pixels in the shadow-covered area of ​​the original remote sensing image sequence based on the terrain shadow directional features to generate a shadow-free remote sensing image includes: Search for the non-shadow reference pixel closest to the pixel in the shadow-covered area in the original remote sensing image sequence along the normal direction of the terrain shadow direction feature; Based on the radiance value of the non-shadow reference pixel, and combined with the spatial distance between the shadow-covered pixel and the non-shadow reference pixel, the reconstructed radiance value of the shadow-covered pixel is calculated by inverse distance interpolation; The reconstructed radiance value replaces the original radiance value of the corresponding pixel in the original remote sensing image sequence to obtain the deshaded remote sensing image.

5. The UAV remote sensing mapping image processing method based on image enhancement according to claim 1, characterized in that, The step of constructing local texture enhancement weights for each pixel based on the gradient magnitude difference between adjacent pixels in the deshaded remote sensing image includes: Calculate the gradient components of each pixel in the deshaded remote sensing image in the horizontal and vertical directions, and obtain the gradient magnitude of each pixel based on the gradient components. Centered on each target pixel, the average gradient magnitude of all neighboring pixels other than the target pixel is calculated within a neighborhood window of a preset radius. The ratio of the gradient magnitude of the target pixel to the average gradient magnitude is used as the local texture enhancement weight of the target pixel.

6. The UAV remote sensing mapping image processing method based on image enhancement according to claim 1, characterized in that, The step of using the local texture enhancement weights to non-linearly stretch the grayscale dynamic range of the de-shaded remote sensing image and outputting the enhanced remote sensing image includes: The normalized gray value is obtained by dividing the original gray value of each pixel in the deshaded remote sensing image by the maximum gray value of the entire image. Using the local texture enhancement weight as the exponent, the normalized gray value is subjected to a power transformation to obtain the enhanced gray value. The enhanced grayscale value is multiplied by the maximum grayscale value of the entire image and then rounded to obtain the output grayscale value of each pixel. The enhanced remote sensing image is then constructed from the output grayscale values.

7. The UAV remote sensing mapping image processing method based on image enhancement according to claim 1, characterized in that, The acquisition of the raw remote sensing image sequence collected by the multispectral sensor on the UAV at multiple heading angles includes: The drone is controlled to fly along a preset circular track within the target area, and at equally spaced angular positions on the circular track, a multispectral sensor is triggered to simultaneously acquire images in the visible light band and the near-infrared band. Visible light and near-infrared images acquired at the same angular location are registered together to form a multi-band original image at that angular location. The original multi-band images at all angular positions are arranged in ascending order of heading angle to form the original remote sensing image sequence.

8. The UAV remote sensing and mapping image processing method based on image enhancement according to claim 2, characterized in that, The step of determining the projection angle of the illumination direction on the horizontal plane at each pixel point based on the sign and magnitude of the rate of change, to form the initial illumination direction field, includes: When the rate of change is positive, the projection angle is assigned the value of the heading angle plus 90 degrees; when the rate of change is negative, the projection angle is assigned the value of the heading angle minus 90 degrees. After normalizing the amplitude to the range of zero to one, multiply it by a preset angle perturbation threshold to obtain the angle perturbation offset. The projection angle is added to the angle perturbation offset to obtain the corrected projection angle, and the corrected projection angles of each pixel are combined to form the initial illumination direction field.

9. The UAV remote sensing and mapping image processing method based on image enhancement according to claim 4, characterized in that, The search for the nearest non-shadow reference pixel to the pixel in the original remote sensing image sequence along the normal direction of the terrain shadow feature includes: Centered on the pixel in the shadow-covered area, two search rays are emitted along the normal direction in the positive and negative directions, respectively. Expand outward pixel by pixel along each search ray until a pixel with a radiance value greater than the preset shadow determination threshold is encountered, and that pixel is taken as a candidate reference pixel. Compare the Euclidean distances between the candidate reference pixels found by the two search rays and the pixels in the shadow-covered area, and select the candidate reference pixel with the smaller Euclidean distance as the non-shadow reference pixel.

10. The UAV remote sensing and mapping image processing method based on image enhancement according to claim 5, characterized in that, After using the ratio of the gradient magnitude of the target pixel to the average gradient magnitude as the local texture enhancement weight of the target pixel, the method further includes: Sort the local texture enhancement weights of each pixel from smallest to largest, and extract the extreme weights that are in the top 5% and bottom 5% after sorting. The extreme value weights are replaced with the minimum and maximum values ​​of the remaining weights to obtain the cropped local texture enhancement weights. The nonlinear stretching is performed by replacing the original local texture enhancement weights with the cropped local texture enhancement weights.

Images