Crop leaf area index estimation method based on unmanned aerial vehicle image processing

Abstract

The invention relates to the technical field of unmanned aerial vehicle remote sensing, in particular to a crop leaf area index estimation method based on unmanned aerial vehicle image processing. The method comprises the following steps: acquiring a remote sensing image of an unmanned aerial vehicle monitoring area; identifying a topographic relief area according to the unmanned aerial vehicle monitoring area remote sensing image, and performing orthoimage reconstruction on the topographic relief area to generate a topographic relief area orthoimage; identifying a terrain shielding area based on an unmanned aerial vehicle monitoring area remote sensing image; image stitching distortion detection is carried out on the terrain occlusion area to obtain image stitching distortion data; performing spectral interpolation on the gap region based on the image splicing distortion data, and determining a plant growth stage of the gap region after the spectral interpolation; and performing lens distortion detection based on the image splicing distortion data to obtain lens distortion data. The precision and accuracy of crop leaf area index (LAI) estimation are improved based on the unmanned aerial vehicle remote sensing technology.

Description

Technical Field

[0001] The present invention relates to the field of unmanned aerial vehicle (UAV) remote sensing technology, and in particular to a crop leaf area index estimation method based on UAV image processing. Background Art

[0002] High-resolution remote sensing equipment equipped with drones can capture high-precision imagery of crop monitoring areas. However, topographical undulation significantly impacts estimation accuracy, especially in complex terrain such as hilly or mountainous areas. Terrain occlusion and distortion can cause spectral distortion or loss in some areas, impacting accurate LAI estimation. Image stitching and distortion processing is inadequate, and distortion during the stitching process is not effectively corrected, leading to stitching errors that further impact the accuracy of the estimation results. Lens distortion and focal length issues are also not effectively corrected. Images with uncorrected lens distortion distort image information, thus affecting LAI estimation. Regarding spectral data processing, traditional methods often lack sophisticated spectral interpolation and data repair measures, making it difficult to accurately recover spectral data missing due to occlusion and distortion, impacting estimation accuracy. Plant growth stage identification also has limitations. Traditional methods often rely on standardized models or simplified spectral signatures, lacking adaptability to different crop types or growth stages, leading to biased estimation results. Processing efficiency is low, often requiring multiple steps of manual intervention and analysis. This is particularly time-consuming for large-scale farmland monitoring or high-frequency monitoring, making it difficult to meet the needs of real-time monitoring and decision-making. Summary of the Invention

[0003] Based on this, it is necessary for the present invention to provide a crop leaf area index estimation method based on drone image processing to solve at least one of the above technical problems.

[0004] To achieve the above purpose, a crop leaf area index estimation method based on UAV image processing includes the following steps:

[0005] Step S1: obtaining a remote sensing image of the area monitored by the drone; identifying the terrain undulating area based on the remote sensing image of the area monitored by the drone, and reconstructing the orthophoto of the terrain undulating area to generate an orthophoto of the terrain undulating area;

[0006] Step S2: identifying terrain-blocked areas based on remote sensing images of the area monitored by the drone; performing image stitching distortion detection on the terrain-blocked areas to obtain image stitching distortion data; performing spectral interpolation on the gap areas based on the image stitching distortion data, and determining the plant growth stage in the gap areas after spectral interpolation;

[0007] Step S3: performing lens distortion detection based on the image stitching distortion data to obtain lens distortion data; identifying focus component damage characteristics based on the lens distortion data; and performing adaptive focal length adjustment based on the focus component damage characteristics to obtain adaptive adjustment focal length data;

[0008] Step S4: performing plant species identification on the orthophoto of the terrain undulating area according to the adaptively adjusted focal length data to obtain plant species data; and estimating the crop leaf area index according to the plant species data and the plant growth stage.

[0009] By acquiring high-resolution remote sensing images and reconstructing orthophotos of areas with undulating terrain, the present invention effectively overcomes the impact of terrain complexity on image accuracy, ensures the geometric accuracy of the images, and further reduces spectral distortion or data loss caused by terrain occlusion or distortion, thereby providing a more reliable basic data layer. By detecting and interpolating the stitching distortion in terrain-occluded areas, not only is the missing data in the original image restored, but the plant growth stage is also accurately determined based on the restored data. This is crucial for estimating the crop leaf area index, improving the accuracy of crop monitoring, avoiding the neglect of occluded and distorted areas in traditional methods, and significantly reducing errors caused by incomplete data. Based on the detection of lens distortion data and adaptive focal length adjustment, not only does it solve the problem of traditional methods failing to effectively correct lens distortion, but it also optimizes image clarity and focus accuracy by dynamically adjusting the focal length, reducing the impact of image distortion. Effective lens distortion correction improves the accuracy of image information, especially in complex terrain areas, and can better display ground details, providing higher-quality data support for subsequent plant species identification and LAI estimation. In terms of plant species identification and growth stage estimation, the improved method overcomes the limitations of traditional methods that rely on simplified models by making full use of restored image data and plant growth information. It can more accurately identify different crop types and their growth stages, reduce deviations caused by unsuitable models, and further improve the accuracy of LAI estimation. By fully processing image stitching and distortion, the processing efficiency in large-scale or high-frequency monitoring is significantly improved, and the need for manual intervention and multiple analyses is reduced, making this method not only suitable for traditional farmland monitoring, but also able to meet the needs of real-time monitoring and decision-making. It effectively solves the impact of factors such as terrain undulation, occlusion, distortion, and lens focal length on the accuracy of LAI estimation, and through precise image restoration, spectral interpolation, plant identification, and focal length adjustment, it significantly improves the accuracy of the estimation results and processing efficiency, providing an efficient and accurate technical path for crop monitoring. BRIEF DESCRIPTION OF THE DRAWINGS

[0010] Other features, objects and advantages of the present invention will become more apparent upon reading the detailed description of non-limiting embodiments thereof made with reference to the following drawings:

[0011] Figure 1 This is a schematic diagram of the steps of the crop leaf area index estimation method based on UAV image processing of the present invention;

[0012] Figure 2 Detailed step flow diagram of step S1 in the present invention;

[0013] Figure 3 Detailed step flow diagram of step S15 in the present invention;

[0014] The purpose, features and advantages of the present invention will be further described with reference to the accompanying drawings and in conjunction with the embodiments. DETAILED DESCRIPTION

[0015] The following is a clear and complete description of the technical method of the present invention in conjunction with the accompanying drawings. Obviously, the embodiments described are part of the embodiments of the present invention, not all of the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without making any creative work are within the scope of protection of the present invention.

[0016] In addition, the accompanying drawings are merely schematic illustrations of the present invention and are not necessarily drawn to scale. Identical reference numerals in the figures denote identical or similar parts, and thus repetitive descriptions thereof will be omitted. Some of the block diagrams shown in the accompanying drawings are functional entities that do not necessarily correspond to physically or logically separate entities. These functional entities may be implemented in software, in one or more hardware modules or integrated circuits, or in different network and / or processor and / or microcontroller approaches.

[0017] It should be understood that although the terms "first," "second," and the like may be used herein to describe various elements, these elements should not be limited by these terms. These terms are used solely to distinguish one element from another. For example, a first element may be referred to as a second element, and similarly, a second element may be referred to as a first element, without departing from the scope of the exemplary embodiments. The term "and / or" as used herein includes any and all combinations of one or more of the listed associated items.

[0018] To achieve this, please refer to Figures 1 to 3 The present invention provides a crop leaf area index estimation method based on UAV image processing, the method comprising the following steps:

[0019] Step S1: obtaining a remote sensing image of the area monitored by the drone; identifying the terrain undulating area based on the remote sensing image of the area monitored by the drone, and reconstructing the orthophoto of the terrain undulating area to generate an orthophoto of the terrain undulating area;

[0020] In this embodiment, a high-resolution remote sensing camera (such as an RGB camera, a near-infrared camera, etc.) carried by an unmanned aerial vehicle is used to take aerial photos of the crop growing area, ensuring that the shooting distance is kept within 2 meters each time, thereby ensuring that the image resolution obtained reaches at least 5 cm / pixel. During shooting, the flight altitude is set to 150 meters to ensure that the coverage area reaches a range of 500 meters × 500 meters. Then, image correction is performed through ground control points (GCPs) to ensure the spatial accuracy of the captured image. Based on the acquired remote sensing image, image processing software (such as Pix4D, Agisoft Metashape, etc.) is used to identify the terrain undulation area. An image segmentation method based on a slope threshold (such as setting the slope threshold to 15°) is usually used to identify the area with larger undulations. For the identified terrain undulation area, an orthophoto reconstruction algorithm (such as structured light, stereo matching algorithm) is used to convert the oblique image into an orthophoto image on the ground that is consistent with the real terrain. During the reconstruction process, the actual coordinates of the ground control points need to be input first to ensure the accuracy of the geographic information of the reconstruction result. The reconstructed orthophoto image should have a spatial resolution of at least 10 cm to ensure a clear display of the terrain features.

[0021] Step S2: identifying terrain-blocked areas based on remote sensing images of the area monitored by the drone; performing image stitching distortion detection on the terrain-blocked areas to obtain image stitching distortion data; performing spectral interpolation on the gap areas based on the image stitching distortion data, and determining the plant growth stage in the gap areas after spectral interpolation;

[0022] In this embodiment, by analyzing the remote sensing images taken by the drone, an edge detection-based algorithm (such as Canny edge detection) is used to identify the terrain obstruction areas in the image. At this time, based on the contrast change of the image, a threshold is set (for example, the edge intensity threshold is set to 50) to determine the area in the image that is obstructed by objects such as trees and buildings. Next, an image splicing distortion detection method (such as reconstruction error detection based on the least squares method) is used to detect the splicing gap area in the image to obtain splicing distortion data. The splicing distortion data can be obtained by calculating the pixel deviation of the overlapping area of the image. Afterwards, based on the splicing distortion data, a spectral interpolation method is used to fill these gap areas. Specifically, a weighted interpolation algorithm, such as the Kriging interpolation method, is used to supplement the missing spectral data by weighted averaging based on the spectral values of adjacent pixels. During the interpolation process, the interpolation window size is set to 5×5 pixels to ensure the accuracy of the calculation. After the spectral interpolation is completed, the vegetation index (such as NDVI, normalized vegetation index) is calculated and combined with the spectral characteristics of the plant to further infer the plant growth stage in each gap area. At this time, an NDVI threshold is set (for example, NDVI>0.3 indicates healthy vegetation, and NDVI<0.3 indicates drought or unhealthy vegetation) to classify the plant growth status of each area and obtain the plant growth stage of the area.

[0023] Step S3: performing lens distortion detection based on the image stitching distortion data to obtain lens distortion data; identifying focus component damage characteristics based on the lens distortion data; and performing adaptive focal length adjustment based on the focus component damage characteristics to obtain adaptive adjustment focal length data;

[0024] In this embodiment, the image stitching distortion data is processed, and the distorted parts in the image are identified using a lens distortion detection algorithm (such as radial distortion detection). By calculating the radial distance from the center of the image to the pixel point, combined with the geometric characteristics of the image, the lens distortion parameters are fitted using the least squares method. The lens distortion data includes distortion coefficients (such as k1, k2, k3, etc.) and principal point offsets, which can be obtained through calibration images. Next, based on the lens distortion data, the image is dedistorted using an image distortion correction algorithm. If focal instability or abnormality is found (for example, the focal adjustment error is greater than 5%), it is identified as a damage to the focusing component. The damage feature can be identified by detecting the degree of focus blur of multiple scenes in the image, and setting a focus blur threshold (for example, a blur greater than 3 pixels indicates an abnormality). After the damage feature is identified, adaptive focus adjustment is performed. According to the current focus setting and image quality, the focus of the drone camera is automatically adjusted to achieve the best imaging effect. During the focus adjustment process, the camera's internal sensors (such as the accelerometer and gyroscope) are used to obtain real-time flight status. Combined with the current flight speed and altitude, the focus range is adjusted to adapt to different shooting distances, and the adjusted focus data is finally obtained.

[0025] Step S4: performing plant species identification on the orthophoto of the terrain undulating area according to the adaptively adjusted focal length data to obtain plant species data; and estimating the crop leaf area index according to the plant species data and the plant growth stage.

[0026] In this embodiment, based on the adaptively adjusted focal length data, plant species identification is first performed on the orthophoto of the terrain undulating area. Using classification algorithms such as support vector machines (SVM) or random forests, different types of plant species are identified by training the spectral features of high-resolution images. During the training process, sample category labels are set (for example, crop A, crop B, weeds, etc.), and species are distinguished by extracting the spectral features of vegetation (such as NDVI, red edge index, etc.). During the recognition process, a 5×5 pixel sliding window is used to extract features to ensure that each plant species area is accurately identified. The recognition result is the plant species category of each pixel. Then, based on the plant species data and the previously calculated plant growth stage, combined with the crop leaf area index (LAI) estimation model, the leaf area index of each species is estimated by polynomial regression or empirical formula. The LAI estimation formula can be performed using the following empirical formula:

[0027]

[0028] Among them, NIR and RED are the spectral values of near infrared and red light bands respectively. max and RED maxis the maximum value, A is the projected area of the crop, and N is the number of pixels. Using this formula, the LAI value of each area is estimated based on species type and growth stage, thus obtaining the overall estimation result of the crop leaf area index.

[0029] Preferably, step S1 is specifically as follows:

[0030] Step S11: Acquire remote sensing images of the drone monitoring area;

[0031] In this embodiment, a rotary-wing UAV equipped with a multispectral imaging module and a high-precision RTK-GNSS module is dispatched to perform low-altitude flight data collection on the target monitoring area. When planning the flight path, grid division is performed according to the boundary coordinates of the target area, and the route overlap is set, wherein the forward overlap is set to 80% and the lateral overlap is set to 70%. The UAV's flight altitude is fixed at 60 meters, the lens's downward angle is set to 90 degrees to obtain vertical remote sensing images, and the shooting interval is set to 1.5 seconds. The flight time needs to be controlled within a time period with sufficient local sunshine and less than 5% cloud cover to ensure the stability of the image spectral information. During the flight, the UAV collects images of five spectral bands: red, green, blue, near-infrared, and short-wave infrared, and the image resolution is set to 0.1 meters / pixel. All images are cached through the UAV's built-in SD card and exported to the ground station computing terminal for further processing after the flight.

[0032] Step S12: generating point cloud data of the monitoring area according to the remote sensing image of the monitoring area by the UAV, and generating a digital elevation model based on the point cloud data of the monitoring area;

[0033] In this embodiment, after the remote sensing image is acquired, the image feature points are extracted using the SIFT (Scale Invariant Feature Transform) algorithm, and the feature point matching relationship between pairs of images is used to generate a three-dimensional sparse point cloud through the Structure from Motion (SfM) method. On this basis, the sparse point cloud is expanded into a high-density point cloud using the Multi-View Stereo (MVS) reconstruction algorithm. The three-dimensional coordinates of the point cloud are solved using GNSS / IMU data and error-aligned with the known ground control point coordinates, with the error controlled within the range of ±5cm. The generated point cloud data format is the LAS standard format, and the point spacing is controlled between 5cm and 10cm. Based on the high-density point cloud data, the DEM (digital elevation model) is reconstructed using the grid interpolation method, and the terrain surface is constructed using the TIN triangulation interpolation method, and the grid is processed into an elevation value matrix with a resolution of 1m×1m, and the output is in GeoTIFF format. The elevation value is measured in meters based on sea level, and the range is set between 0m and 200m to ensure coverage of all terrain features in the monitoring area.

[0034] Step S13: extracting regional pixel height values according to the digital elevation model, and calculating the slope of the monitoring area based on the regional pixel height values; estimating spectral reflectance values according to the digital elevation model, and identifying shadow areas based on the spectral reflectance values;

[0035] In this embodiment, based on the elevation value of each grid point in the digital elevation model, a five-point difference method is used to calculate the difference between the adjacent elevation values of each pixel in the east-west and north-south directions, respectively, to obtain the slope change trend in these two directions. By comparing the elevation differences of adjacent grid points and combining them with the actual resolution, the ground slope angle value of each pixel is calculated. All slope angle values are output in a grid format to form a slope map. To ensure the accuracy of the slope calculation, a grid window with a fixed step size of 1 meter is used for processing. The difference window extends two pixels in all directions in the elevation matrix and slides pixel by pixel. In the resulting slope map, areas with slope angles greater than 15 degrees are defined as high slope areas, and the slope level is labeled as an integer ranging from 0 to 90 degrees. Areas with slopes less than 3 degrees are classified as flat areas. This classification result serves as the basis for subsequent analysis of undulating terrain. In the shadow area identification process, the acquired multispectral remote sensing image is first atmospherically corrected, and the reflectivity of each spectral band in the remote sensing image is corrected using the 6S radiation transfer model. Required parameters include the measured solar altitude angle of 45 degrees at the time the remote sensing image was captured; the measured average ground reflectance of 0.25; a continental aerosol model with an aerosol optical depth of 0.1. All parameters are derived from synchronized data from on-site optical environmental monitoring equipment. After calibration, the reflectance values for each pixel in the red and near-infrared bands are extracted, and the vegetation index value is calculated through band interpolation and normalization. Pixels with an NDVI less than 0.05 and a total radiation intensity in the visible band less than 300 watts per square meter are marked as shadow pixels. This intensity value is measured by a ground-based solar radiation sensor and compared with the spectral intensity of the remote sensing image, setting a static threshold. Shadow identification results are output as a shadow mask image in GeoTIFF format, where a pixel value of 1 indicates shadow and a value of 0 indicates non-shadow. This image is named "shadow_mask.tif" and will be used for spatial overlay analysis in the subsequent terrain relief identification step.

[0036] Step S14: determining the terrain undulation area based on the slope of the monitoring area and the shadow area;

[0037] In this embodiment, the slope map and the shadow area mask map are combined to identify the terrain undulating area through a dual discrimination strategy. First, the area with a slope value greater than 15° is set as the first type of undulating area candidate mask (mask_slope), and the shadow area identification mask is set as the second type of candidate mask (mask_shadow). Then, a logical "and" operation is used to generate an intersection mask (mask_combined), that is, only when a pixel meets the conditions of a slope greater than 15° and being marked as a shadow, the pixel is classified as undulating terrain. In order to exclude non-continuous small area noise points, morphological dilation processing (kernel is 3×3) and 8-connected domain extraction are performed on the intersection mask to filter out isolated areas with less than 200 pixels. Finally, a binary mask of the terrain undulating area (elevation_mask.tif) is formed, and the boundary of its covered area is converted by geographic coordinate projection and recorded in vector format (shapefile) for subsequent ortho reconstruction.

[0038] Step S15: reconstructing the orthophoto of the terrain undulating area to generate an orthophoto of the terrain undulating area.

[0039] In this embodiment, the mask map of the terrain undulation area is used as input, and the original remote sensing image is cropped to obtain the ROI area image set. Based on these images, the Bundle Adjustment technology based on multiple images is used to jointly solve the internal parameters (focal length, principal point coordinates) and external parameters (camera attitude and position) between the images, and the GNSS / IMU and ground control point data are strictly aligned, and the error is controlled within 2 pixels. During the reconstruction process, the RFM (Rational Function Model)-based geometric model is used for image correction to eliminate the tilt distortion and parallax in the image. The resampling resolution of the corrected image is set to 0.1 meters, and the unified projection is the WGS-84UTM coordinate system. The final synthesized orthophoto output format is GeoTIFF, and the file naming rule is "orthophoto_hilly_zone_date.tif". The sensor parameters, shooting time, route number and area code are recorded in the image metadata as input for subsequent vegetation recognition and LAI estimation.

[0040] Preferably, step S15 is specifically as follows:

[0041] Step S151: collecting remote sensing images of the terrain undulating area to obtain remote sensing images of the terrain undulating area;

[0042] In this embodiment, a multi-rotor vertical take-off and landing drone equipped with a multispectral imager is used to capture remote sensing images of areas with undulating terrain. The route design requires covering the entire boundary of the target area with undulating terrain and overlapping with a redundant buffer of at least 20 meters to avoid missing edge data. The flight altitude is set to 80 meters, the heading overlap rate is set to 85%, and the lateral overlap rate is set to 75%. The flight time is selected to complete the shooting within half an hour before and after local noon to avoid long shadows interfering with the image quality. The imaging equipment requires a shooting resolution of 0.05 meters / pixel and a shooting time interval of no more than 2 seconds. Multispectral images include five bands: red, green, blue, red edge, and near-infrared. After the original image is collected, it is stored in .tiff format, and the image number is automatically encoded according to the flight order.

[0043] Step S152: Acquire ground control point data; extract corner points of the remote sensing image of the terrain relief area, wherein the number of extracted corner points is set to 50-100; perform local feature point registration based on the ground control point data and the corner points to obtain local registration data;

[0044] In this embodiment, the ground control points are laid out and data collected using RTK differential GPS measurement equipment. The number of control points is no less than 8, distributed around the target area and the central area. The longitude, latitude and altitude values are recorded for each control point, and the control point error is controlled within ±0.05 meters. The corner point extraction operation is performed using the Harris corner point detection method, with the corner point response threshold set to 0.01 and the non-maximum suppression window set to 3×3. The number of corner points extracted in each remote sensing image is between 50 and 100, and the specific number is determined based on the complexity of the image texture. After the corner point extraction is completed, the corner points are feature matched based on the geographic coordinates of the control points, and the RANSAC algorithm is used to eliminate mismatched points, and the matching residual is controlled within ±0.3 pixels. The final output local registration data includes a one-to-one correspondence between the image coordinates and the actual coordinates of the corner points in the image, which is saved in *.txt format for subsequent image registration conversion operations.

[0045] Step S153: performing image projection conversion on the orthophoto of the terrain relief area according to the local registration data to obtain a local orthophoto projection;

[0046] In this embodiment, based on the local registration data, bilinear interpolation is used to reconstruct the spatial position of the image pixels and perform projection transformation. The original coordinate system of the image is the image pixel space (Image Space), the target projection coordinate system is the WGS84 / UTM projection system, the area code is determined according to the longitude where the drone image is taken, and the resolution remains unchanged at 0.05 meters / pixel. During the projection transformation process, coordinate reprojection processing is performed on each image separately, and the conversion accuracy control threshold is set to 0.001 meters during the processing. The output result is a local orthographic projection image. The image coordinate system has complete spatial reference information. The file format is GeoTIFF, and the spatial projection parameters are embedded in the header file.

[0047] Step S154: performing texture correction on the local orthographic projection to obtain a texture-corrected local orthographic projection;

[0048] In this embodiment, the color and texture of the local orthophoto projection image are uniformly adjusted. First, the image histogram distribution is calculated, and the grayscale histogram mean and standard deviation of each image are extracted. The image grayscale is adjusted using histogram matching technology. The target image is the image with the grayscale mean closest to the median in the entire image sequence, and the standard deviation is controlled within ±5 grayscale levels. To eliminate the impact of illumination differences on texture consistency, the linear illumination normalization method is further used to adjust the brightness values of the red, green, and blue bands. After adjustment, the brightness standard deviation of the image in the three bands is controlled within ±3 grayscale levels. After processing is completed, the texture-corrected local orthophoto is output and named "corrected_tile_number.tif".

[0049] Step S155: performing image stitching based on the texture-corrected local orthophoto projection to generate an orthophoto of the terrain relief area.

[0050] In this embodiment, before performing the stitching process, the edge overlapping areas of all texture-corrected local orthophoto images are first extracted, and feature point extraction and SIFT matching operations are performed based on the overlapping pixel areas in the middle of the images. At least 120 sets of feature point pairs are matched for each pair of images, and the weighted average method is used to calculate the brightness and color transition boundaries of the overlapping areas. The stitching algorithm adopts an image registration and fusion method, and uses a weighted fusion method to perform brightness transition processing on the overlapping areas. The fusion coefficient of overlapping pixels is set to 0.6:0.4 to ensure that the transition of the stitched image boundaries is smooth and there are no obvious seams. The image output is stitched using a single coordinate system and consistent resolution, and finally a unified orthophoto image of the terrain undulating area is generated. The output format is GeoTIFF, the file is named "terrain_ortho.tif", and it is embedded with complete spatial projection information and image stitching metadata.

[0051] Preferably, the identification of the terrain obstruction area in step S2 is specifically as follows:

[0052] Grayscale conversion is performed based on the remote sensing image of the UAV monitoring area to obtain a grayscale remote sensing image;

[0053] In this embodiment, when grayscale conversion is performed on remote sensing images of the drone-monitored area, the grayscale values are calculated using a weighted average of the red, green, and blue bands in the remote sensing image. The weights for the red, green, and blue bands in the grayscale conversion formula are set to 0.299, 0.587, and 0.114, respectively, to ensure that the converted grayscale image accurately reflects the intensity variation characteristics of the target area. The image pixel resolution is uniformly set to 0.05 meters per pixel, and the image size is maintained within the original remote sensing image acquisition size range. The output image format is an uncompressed grayscale 8-bit TIFF format. No image sharpening, edge enhancement, or contrast stretching is performed during this process to preserve the original texture characteristics.

[0054] Calculate the gray-level co-occurrence matrix of gray-level remote sensing images;

[0055] In this embodiment, after the grayscale image is generated, the grayscale co-occurrence matrix of the grayscale remote sensing image is calculated, and four main directions (0°, 45°, 90°, 135°) are selected, and the number of grayscale levels is set to 64. The window size is set to a sliding window of 21×21 pixels with a step size of 10 pixels to ensure that local texture features can be stably extracted. The statistical distance parameter in the grayscale co-occurrence matrix is set to 1 pixel, and the calculation method is to perform joint probability statistics on the pixel grayscale pairs in each sliding window, construct the co-occurrence matrix in the corresponding direction, and perform normalization. Each window generates a set of grayscale co-occurrence matrices in four directions, which are stored in separate data structures to support the calculation of subsequent statistical features.

[0056] According to the gray-level co-occurrence matrix, the roughness is calculated; According to the gray-level co-occurrence matrix, the entropy value is calculated;

[0057] In this embodiment, after completing the construction of the grayscale co-occurrence matrix, the roughness and entropy of the image texture are respectively counted according to the above co-occurrence matrix. The roughness calculation formula adopts a combination of variance and contrast indicators. The variance is obtained by accumulating the square difference of the grayscale value from the mean, and the contrast is achieved by weighted summation of the square difference of the diagonal elements of the co-occurrence matrix. The entropy value calculation adopts the information entropy definition formula. The normalized value in each co-occurrence matrix is -log transformed and multiplied by itself and summed. The roughness threshold is set to 0.12, and the entropy threshold is set to 2.5. The two thresholds are based on the statistics of the field sample image of the target area.

[0058] Identify uneven texture areas in grayscale remote sensing images based on roughness and entropy values;

[0059] In this example, if the roughness value of a sliding window is greater than 0.12 and the entropy value is greater than 2.5, the area corresponding to the window is marked as an uneven texture area. This operation generates a label map in the image dimension, where each pixel position stores whether it belongs to an uneven texture area. The output format is a Boolean mask map (1 indicates uneven texture, 0 indicates uniform texture area).

[0060] Detect linear transition edges based on remote sensing images of the drone monitoring area; identify ridge areas in uneven texture areas based on linear transition edges; obtain the sunlight angle; simulate the illumination of the ridge area based on the sunlight angle, and identify terrain occlusion areas during the simulation process.

[0061] In this embodiment, after the uneven texture region is identified, the Canny edge detection method is used to identify edges in the remote sensing image of the drone monitoring area. The high and low thresholds in the Canny algorithm are set to 0.05 and 0.15, respectively, and the 3×3 Sobel operator is used to calculate the gradient amplitude and direction. To identify linear transition edges, the edge map is subjected to a Hough transform, with a detection line segment length threshold of 30 pixels and an angular resolution of 1 degree. Edge segments with a linear fitting error of less than 2 pixels are defined as linear transition edges. The detected linear edges are superimposed and analyzed with the uneven texture region mask. If a linear edge lies on the boundary of the uneven texture region and its direction is consistent with the local texture direction (the angle difference is less than 10 degrees), the edge is marked as a ridge region line. To obtain the solar illumination angle, the specific time information of the image time point is collected, and the solar altitude and azimuth angle data of the region are queried. The calculation is based on a public solar radiation angle database combined with the image capture time and geographic coordinates. Taking June 10, 2024, at 12:15 noon at 23.5°N and 113.2°E as an example, the solar altitude angle is 75° and the solar azimuth angle is 145°. This solar angle parameter is input into the digital elevation model, and the terrain is simulated using the illumination simulation algorithm (Hillshade algorithm). The light source direction is incident according to the solar angle. The orientation and slope of each terrain pixel are calculated and compared with the solar incidence direction. If the angle between the incidence direction and the slope direction exceeds 90°, the pixel is blocked. This operation outputs a mask map of the terrain blocked area. After superimposing it with the ridge area map, the final annotation map of the light blocked area is obtained, which is used to eliminate the influence of unreasonable illumination areas in the subsequent leaf area index calculation.

[0062] Preferably, the identification of the terrain obstruction area in step S2 is specifically as follows:

[0063] Collect multiple remote sensing images of terrain-blocked areas to obtain multiple remote sensing images;

[0064] In this embodiment, in the process of acquiring multiple remote sensing images of the terrain-obstructed area, a drone is used to repeatedly fly over the target area at different times and angles to collect image data. The drone's flight altitude is set to 80 meters, the shooting overlap rate is controlled to 80% in the forward direction and 70% in the lateral direction, the camera is fixed in a ground-facing vertical downward shooting mode, and it is ensured that each flight mission records complete GPS coordinate information and shooting timestamps. No less than 3 sets of image sequences are collected each time, the time interval is controlled within 2 hours, and the illumination change does not exceed ±15 degrees to ensure stable image illumination conditions. The image format is required to be RAW format with a resolution of 5472×3648 pixels to ensure that the detail resolution can support subsequent occlusion structure analysis.

[0065] Generate a terrain occlusion splicing image based on multiple remote sensing images; identify terrain occlusion edge segments based on the terrain occlusion splicing image; calculate the line segment slope of the terrain occlusion edge segment; identify the angle jump line segment of the terrain occlusion edge segment based on the line segment slope;

[0066] In this embodiment, based on the multiple remote sensing images acquired as described above, the image is spliced using a feature point matching method to generate a terrain occlusion spliced image. The SIFT feature extraction algorithm is used to extract no less than 5,000 feature points from each image, and the FLANN matching method is used to pair the feature points between images. Feature point pairs with a matching error of less than 1.5 pixels are selected for homography matrix estimation. The RANSAC algorithm is used to eliminate mismatched point pairs, and the threshold for the homography matrix fitting error is set to within 2 pixels. Multiple images are spliced using a geometric registration method, and a multi-band fusion method is used to achieve smooth transitions in the image boundary areas. The final output resolution remains unchanged from the original image. In the spliced terrain occlusion image, edge segments are extracted using the Canny edge detection method, and edge segments are detected using the Hough line detection method. The minimum length of the segment is set to 20 pixels, and the angle accuracy is 1 degree. A linear slope is fitted to each edge segment, and the slope calculation is performed using the difference division of the coordinates of the two end points (k = Δy / Δx). Segments with Δx < 5 pixels are filtered to avoid excessive slopes that affect the analysis. Analyze the slope changes between consecutive line segments. If the absolute value of the slope difference between two adjacent line segments is greater than 1.2 and the line segment angle changes by more than 25 degrees, it is determined to be an angle jump line segment, and its start and end pixel coordinates are recorded as the abnormal line segment point set.

[0067] Extract terrain occlusion connected areas based on terrain occlusion stitching images; search for 8 adjacent edge pixels based on terrain occlusion connected areas, and perform edge chain extension based on the 8 adjacent edge pixels to obtain edge chain data; detect edge chain breakpoints based on the edge chain data;

[0068] In this embodiment, when extracting terrain-occluded connected areas in a spliced image, a Flood Fill algorithm is used to extract connected areas based on a binary occlusion mask map. The pixel adjacency method adopts an 8-adjacency mode, and the minimum area of a single connected area is set to 500 pixels. All edge pixel coordinates are extracted on the boundary of each connected area, and each pixel is queried for whether there is an occlusion mask edge pixel among its 8 adjacent pixels. If so, the edge chain is continued to be expanded. The edge chain expansion process adopts a depth-first traversal method until no adjacent pixels can be expanded. All edge chain results are saved in the form of a coordinate array. Breakpoint detection is performed on each edge chain by calculating the Euclidean distance between consecutive points in the edge chain. If the distance is greater than 3 pixels and there are no other adjacent points connected, it is marked as a breakpoint, and its index position and adjacent coordinates are recorded.

[0069] Determine the geometric structure distortion area image of the terrain occlusion stitching image based on the angle jump line segment and the edge chain break point;

[0070] In this embodiment, the aforementioned angle-jump segment data and edge chain breakpoint data are combined for spatial positional matching. If the breakpoint is less than 5 pixels away from the angle-jump segment endpoint, or if they are located in the same region (with the same connected domain number), they are merged into a geometrically distorted region. The region coordinate bounding box is output and stored as a mask. This geometrically distorted image is numbered to identify different regions, and this boundary information is used for subsequent posture data mapping analysis.

[0071] Extract image acquisition timestamp based on terrain occlusion stitching image; obtain UAV three-axis angle data; map image acquisition timestamp to UAV three-axis angle data to obtain image acquisition UAV three-axis angle data;

[0072] In this embodiment, the image acquisition timestamp is extracted from the stitched image. The field "DateTimeOriginal" is obtained by reading the original metadata (EXIF) of the image, and the extracted time format is "YYYY:MM:DD HH:MM:SS". The data record row that matches the timestamp is searched in the drone sensor record file. The file contains three-axis angles (pitch angle Pitch, roll angle Roll, yaw angle Yaw), and a timestamp index table is established with a time accuracy of 0.01 seconds. After matching, the corresponding three-axis angle data is extracted and stored as a triplet, for example (Pitch = 2.5°, Roll = -0.8°, Yaw = 175.3°), and the unit is angle.

[0073] Based on the image acquisition of the UAV's three-axis angle data, the UAV attitude jump detection is performed to obtain the UAV attitude jump data;

[0074] In this embodiment, the attitude change rate is calculated in time series from the three-axis angle data corresponding to all image acquisition timestamps. The calculation formula is the difference between the three-axis angles of the current frame and the previous frame divided by the time interval. If the pitch or roll change rate of a frame exceeds 15° / second, or the yaw change rate exceeds 25° / second, it is determined to be a drone attitude jump frame, and the image frame number and acquisition time are recorded.

[0075] Extracting attitude jump image frames based on the attitude jump data of the UAV; mapping the attitude jump image frames to the terrain occlusion stitching image to obtain the attitude jump terrain occlusion image;

[0076] In this example, the corresponding image frames are extracted and aligned with the spatial locations in the terrain-occluded mosaic image. Based on the correspondence between their GPS coordinates and the coordinates in the mosaic image, they are mapped to the mosaic image. This mapping method uses an affine transformation approach used in image registration, performing a linear transformation between image corner points and spatial points in the mosaic image to produce a posture-jumping terrain-occluded image.

[0077] The posture jump terrain occlusion image and the geometric structure distortion area image are integrated to obtain the image stitching distortion data.

[0078] In this embodiment, all attitude-jump terrain occlusion images and geometric distortion region images are spatially overlapped. The two image masks are overlaid at the pixel level, and the pixel positions of the intersection area are counted as image stitching distortion data. The system also outputs information such as region number, start and end coordinates, and distortion type flags (angle jump, chain break, attitude jump). The results are stored as multi-channel image masks and structured table data in CSV format, which are used to eliminate abnormal image areas when calculating the leaf area index.

[0079] Preferably, the slit region spectrum interpolation in step S2 is specifically as follows:

[0080] Extract zero pixel data based on image stitching distortion data; locate gap areas based on zero pixel data;

[0081] In this embodiment, after the image is acquired, the multispectral image (including red, green, blue, red edge, and near-infrared bands) taken and spliced by the drone and the corresponding image stitching mask data are first read. The image matrix reading function in OpenCV is used to traverse the entire image pixel by pixel to determine whether the pixel values of all channels in the image are 0 at the same time. If they are 0, they are considered invalid pixels caused by stitching distortion. The two-dimensional coordinate information of these zero-value pixels is recorded as a coordinate list. To ensure processing efficiency, only areas with a zero pixel ratio of no more than 5% in the image are processed. The excess area needs to reprocess the image stitching parameters before proceeding with this process. After obtaining all zero-pixel coordinate points, the eight-adjacent region connection method is used to aggregate them. The connected region labeling method in the morphological operation is used, and the scipy.ndimage.label function is used in the Python environment to cluster adjacent zero pixel points into connected regions. The minimum pixel number threshold of the connected region is set to 25, and areas below this value are regarded as isolated points and ignored. The bounding information of the circumscribed rectangle of the eligible area is extracted, and the position index of each gap area is recorded for subsequent interpolation positioning.

[0082] A 5x5 neighborhood window is used to extract the effective spectrum value of the gap area;

[0083] In this embodiment, the pixel coordinate point of each gap area is used as the center and expanded to the surrounding area to form a 5×5 neighborhood sliding window with a fixed side length of 5 pixels. The image boundary filling technology is used to process the edge area to ensure that the neighborhood can be successfully extracted for each pixel. For the remaining 24 pixels in the window except the central pixel, all channel spectral values in the multispectral image are read respectively. The validity standard is that the spectral value is non-zero and the reflectivity is between 0.05–0.95. The neighborhood with less than 15 valid pixels does not participate in the interpolation process and is recorded as an unrepairable area. The valid pixels and their distance to the center point are retained for weighted processing.

[0084] Performing weighted average calculation based on the effective spectrum value to obtain a weighted spectrum value; performing spectrum interpolation on the gap area based on the weighted spectrum value to obtain a spectrum interpolation gap area;

[0085] In this embodiment, for each gap pixel to be repaired, the spectral reflectance and corresponding distance weight of the valid pixels in its neighborhood are read, and the weight is calculated using the inverse square of the distance. A weighted average operation is performed on each band, and the output result is the interpolated spectral value of the gap pixel. This operation is implemented in Python using NumPy arrays for matrix processing, improving batch interpolation efficiency. This operation is repeated for all gap pixels in the image, and finally a spectral interpolation map with the same size as the original image is formed. All interpolation areas are limited to zero-pixel areas, and the rest of the original image remains unchanged.

[0086] Calculating vegetation index based on spectral interpolation gap area; identifying vegetation area in spectral interpolation gap area according to vegetation index;

[0087] In this example, the vegetation index is calculated pixel by pixel using the spectrally interpolated image data, selecting the red and near-infrared reflectance data as input. The resulting image is used to generate an NDVI image (normalized vegetation index map). The criterion for identifying vegetation areas is set to an NDVI value greater than or equal to 0.3. All pixels that meet this criterion are marked, and a binary mask image of the vegetation area is generated. This mask image will be used for subsequent leaf spectral sampling and analysis.

[0088] Collect leaf spectra of vegetation areas and invert chlorophyll content based on the leaf spectra;

[0089] In this embodiment, 10 core pixels with NDVI values higher than 0.6 are selected from the vegetation area as representative areas. At the position of each core pixel, the reflectance values of all bands in the original multispectral image are extracted. Each pixel is regarded as a spectral sample, forming five-dimensional vector data containing five bands: red, green, blue, red edge and near infrared. The spatial position and band values of the sampling points are uniformly saved in the sample data structure for use in the inversion chlorophyll analysis stage. Using the sampled leaf spectral values, the vegetation spectral index of each sample is calculated according to the conventional two-band combination method. The selected vegetation spectral index is the ratio of TCARI to OSAVI. The numerical range of the ratio is pre-set to correspond to the measured chlorophyll concentration range, and the numerical mapping is performed by a lookup table. The numerical range is set as follows: a ratio of 0.1–1.2 corresponds to a chlorophyll concentration of 5–45 micrograms per square centimeter, and the mapping relationship is calculated using linear interpolation. The results of all samples are saved in a list format, and their spatial coordinate information is recorded for regional analysis.

[0090] Determine plant growth stage based on chlorophyll content.

[0091] In this embodiment, the chlorophyll concentration data is divided into four growth stages according to the numerical range: less than 15 micrograms per square centimeter is defined as the seedling stage, 15 to 30 micrograms is defined as the vegetative growth period, 30 to 40 micrograms is defined as the jointing and booting period, and more than 40 is the filling and maturity period. The image partitioning statistics method is used to calculate the number and proportion of chlorophyll samples in each stage in the entire image. Combined with the vegetation area distribution mask, the dominant growth stage is determined for each area block. The final output results include: a spatial distribution map of the vegetation area, a growth stage classification map, and an area ratio table corresponding to each stage. This result serves as the input basis for the subsequent leaf area index estimation model.

[0092] Preferably, step S3 is specifically as follows:

[0093] Step S31: performing geometric distortion identification based on the image stitching distortion data to obtain geometric distortion data;

[0094] In this embodiment, during the image stitching process, the control point residual error information and the image edge stretch ratio are obtained. By analyzing the spatial geometric position changes of the feature points before and after stitching, the geometric distortion features present in the image are identified. The SIFT algorithm is used for feature point extraction, and the control point matching uses model inlier screening based on RANSAC (random sampling consistency) to ensure the stability of the error value calculation. The geometric distortion judgment threshold is set to 0.8 pixel error mean. That is, when more than 80% of the feature control points have an offset greater than 0.8 pixels after stitching, it is marked as having geometric distortion. The radial stretch ratio of the image edge to the center point is further used for comparative analysis to calculate the radial scale deviation. When the scale difference in a certain direction is greater than 5% (that is, the difference between the edge extension and the center pixel spacing is greater than 5%), the presence of geometric stretch distortion in that direction is recorded. Finally, the image coordinate range of the geometric distortion area and the error value list are output to form a structured geometric distortion data set.

[0095] Step S32: performing lens damage detection based on the image stitching distortion data to obtain lens damage data;

[0096] In this embodiment, a differential analysis is performed between the light spot intensity distribution diagram in the stitched image and the uniform illumination test diagram. First, before the drone collects images, a standard uniform illumination plate is installed on it for illumination shooting, and a uniform brightness image is obtained as a reference template. After the actual captured image is grayscaled, the pixel difference is calculated with the reference image to generate a brightness difference map. The lens abnormality recognition threshold is set to a brightness difference greater than 15% (that is, the grayscale value of a certain area in the actual image differs from the grayscale of the pixel at the same position in the reference image by more than 15%). When a bright or dark spot area with a mainly circular shape and blurred edges appears in the image and lasts for more than 100 pixels, combined with the spatial position of the area in the center or periphery of the lens, it is judged that the lens has surface damage or coating peeling. The lens light loss area is detected for images of different bands, and the coordinate information of the lens damage area and its area and brightness difference percentage are output to form lens damage data.

[0097] Step S33: Integrate the geometric distortion data and the lens damage data to obtain lens distortion data;

[0098] In this embodiment, the geometric distortion data output in step S31 and the lens damage data output in step S32 are structured and integrated. A unified data format is defined, including the following fields: distortion type (geometric distortion, brightness distortion), coordinates of the center point of the distortion area, area of the distortion area, distortion direction (if any), and distortion intensity value (expressed as pixel error or brightness deviation). A coordinate matching method is used to perform spatial intersection analysis on the two types of distortion data to determine whether there are collinear or overlapping areas. If the center distance between the two types of distortion areas is less than 20 pixels and the overlapping area exceeds 50%, it is marked as a composite distortion area, and the priority of the area is set to high-risk lens distortion. The final output lens distortion data structure contains all independent and composite distortion area information, which serves as an important basic input for subsequent analysis of focal length component damage.

[0099] Step S34: identifying focus assembly damage features based on the lens distortion data;

[0100] In this embodiment, cluster analysis is performed on the composite distortion area in the lens distortion data, and DBSCAN (density-based spatial clustering algorithm) is used to identify multiple high-density distortion blocks, and their spatial distribution and change trends are analyzed. The minimum number of cluster samples is set to 5, and the neighborhood radius is 30 pixels. If the cluster centers of the composite distortion area are distributed on the same side of the image or the offset in a certain direction is too large (the offset angle exceeds 15 degrees), it is judged that the optical axis in the lens is tilted, which is a sign of imbalance in the focus component. Further, by comparing the changes in the degree of distortion in images at different shooting angles, it is determined whether the focus motor has a slow response or position offset. The drift rate of the distortion center in multiple consecutive images is calculated. When the drift exceeds 10 pixels and continues to appear for more than 3 frames, it is recorded as a suspected focus motor failure. All identified damage features are encoded and output in the form of damage position (coordinates), drift value, change trend, etc., providing a basis for focal length adjustment.

[0101] It is particularly important that step S34 includes the following steps:

[0102] Step S341: performing distortion type classification based on the lens distortion data to obtain lens barrel distortion data and lens pincushion distortion data;

[0103] In this embodiment, image data captured by a drone equipped with high-resolution remote sensing equipment is required. This image data contains the effects of lens distortion. Through geometric correction and distortion analysis of the image, image processing algorithms, such as Zhang's calibration method, are used to perform distortion correction analysis on the image. Specifically, the characteristics of lens distortion are first extracted. The image distortion parameters (including radial distortion coefficients and tangential distortion coefficients) are then obtained using a checkerboard calibration method or other suitable calibration methods. These parameters are then used to correct the image distortion. With barrel distortion, the image edges are stretched outward, while with pincushion distortion, the image edges are contracted inward. Using the characteristics of these two types of distortion, the image data can be divided into two categories: barrel distortion data and pincushion distortion data. This division is based on the edge curvature of the distorted image. A threshold is set to determine the distortion characteristics of each pixel. Barrel distortion and pincushion distortion can be determined by calculating the curvature of the image edge. A curvature threshold (e.g., 0.01) is typically used for classification. The curvature of the edge of a barrel-distorted image is greater than a threshold, while the curvature of the edge of a pincushion-distorted image is less than a threshold. Based on this threshold, the image data is divided into two categories: barrel-distorted data and pincushion-distorted data. The subsequent steps will process these two categories of data separately.

[0104] Step S342: identifying focus lag characteristics based on lens barrel distortion data;

[0105] In this embodiment, the barrel-distorted image data obtained in step S341 is extracted and used to identify focus lag. Focus lag is primarily manifested as image blur and inability to accurately focus in certain areas. For barrel-distorted images, an edge detection algorithm, such as Canny edge detection, is used to determine the in-focus areas within the image. Because barrel distortion typically causes image edges to appear extended, focus lag is determined by calculating the change in sharpness across various regions within the image. Specifically, a sharpness evaluation algorithm (such as a Laplace transform) is used to analyze the image for clarity, and a clarity threshold is set (for example, areas with a clarity value less than 5 are considered blurry). Next, areas within the image whose clarity values do not meet the standard are identified; these areas are blurry due to focus lag. By comparing the image with a known standard focus image, the areas caused by focus lag can be more accurately determined. After identifying these areas, the corresponding coordinate information and blur level are recorded to provide a basis for subsequent adaptive focus adjustment.

[0106] Step S343: identifying focus blur characteristics based on the lens pincushion distortion data;

[0107] In this embodiment, the pincushion distortion image data obtained in step S341 is extracted and focus blur is identified. Images with pincushion distortion typically exhibit a sharp center region and blurred edges. The edge distortion prevents accurate focus. To identify focus blur characteristics, frequency domain image analysis methods, such as Fourier transform, can be used to convert the image from the spatial domain to the frequency domain. The spectrum of blurred areas typically contains low-frequency information and has small amplitude values. By calculating the amplitude and phase spectrum of the image in the frequency domain and setting a threshold (for example, frequency components with spectrum amplitudes below a certain threshold, such as 0.1), focus blur areas can be identified. For these areas, a clarity analysis algorithm (such as a root mean square error or gradient method) can be used to further confirm whether focus blur is present. Based on information such as the size and location of the blurred area, its behavior under the influence of lens pincushion distortion is further analyzed, ultimately confirming the specific location and degree of focus blur. This method effectively demarcates focus blur areas and identifies focus blur characteristics caused by lens pincushion distortion.

[0108] Step S344: Integrate the focus lag feature and the focus blur feature to obtain the focus component damage feature.

[0109] In this embodiment, for each identified focus lag area and focus blur area, the spatial overlap and distance calculation methods are used to determine whether these areas have commonalities or influence each other. Specifically, based on the spatial position relationship, the coordinates of the focus lag area and the focus blur area are merged, and the adjacent area analysis is used to detect whether there are repeated blurred areas caused by the failure of the focus component. By comparing the focus clarity and blur changes in different areas, it is determined whether the focus component is damaged. Statistical analysis methods such as weighted averaging or logistic regression are used, combined with the focus lag features and focus blur features, to evaluate the degree of influence of the focus component damage. If the focus lag area and the focus blur area are highly overlapped in space and the clarity difference is large, it can be determined that the focus component is damaged. Finally, these damage features are extracted to provide reference data for the subsequent adaptive focal length adjustment step.

[0110] Step S35: performing adaptive focal length adjustment based on the damage characteristics of the focusing component to obtain adaptively adjusted focal length data.

[0111] In this embodiment, after identifying the specific damage characteristics of the focus component, focus compensation is performed by controlling the autofocus module of the drone gimbal system. The focus adjustment instruction parameters include the adjustment direction (pull in or push out), the adjustment amplitude (in microns, with an accuracy of 10μm), and the adjustment frequency (the default interval is 1 time per second). According to the optical axis offset direction and the distortion drift trend in the damage characteristics, the correction vector of the point where the compensation direction coincides with the target focus is set. For example, if it is found that the upper part of the image continues to drift in distortion, the focus is automatically compensated upward and the focal plane is pulled away. After each adjustment, the image is recaptured and the degree of distortion is re-detected. When the center displacement of the distortion is stable within 5 pixels and the distortion brightness deviation is less than 10%, the current focal length value is recorded as the adaptive focal length in this scene. Finally, the adaptive adjustment focal length data containing the focal length values before and after adjustment, adjustment parameters, number of corrections, etc. is formed and written into the drone parameter log for archiving.

[0112] Preferably, step S31 is specifically as follows:

[0113] Step S311: extracting image corner points based on the image stitching distortion data, and performing corner point matching based on the image corner points to obtain matching corner point data;

[0114] In this embodiment, image corners are extracted from the stitched, distorted data of the stitched images. First, grayscale normalization is performed on each image, and pixel values are uniformly mapped to the range 0–255. Edge smoothing is performed using a Gaussian filter with a kernel size of 5×5 and a standard deviation of 1.2 to reduce the impact of image noise. Corner detection uses the Harris corner detection algorithm, with a response function threshold of R>100,000 as the valid corner determination criterion. The image block window size is 3×3, and the corner sensitivity parameter k is 0.04. After detection, no fewer than 300 valid corners are extracted from each image, and their pixel coordinates and response values are output. In the corner matching stage, a mutual information matching method based on local image blocks is used. For each corner point, an 8×8 grayscale subwindow is extracted as a descriptor. The mutual information score in the matching image is calculated, and the top three candidate points with the highest scores are selected. A secondary screening is performed using a combined score of Euclidean distance and gradient direction angle similarity to select the final matching point pairs. In order to eliminate mismatched data, a geometric consistency constraint based on RANSAC was introduced, the number of iterations was set to 200, the inlier error threshold was set to 1.5 pixels, and finally more than 150 groups of stable matching corner point pairs were retained as matching corner point data.

[0115] Step S312: identifying a transformation matrix relationship based on the matching corner point data;

[0116] In this embodiment, the transformation relationship between images is identified by matching corner point data, and the transformation model is set as the homography matrix (Homography Matrix), and the least squares method is used for solution. Each set of matching corner point coordinate pairs is set as input, and a linear transformation equation group Ax=b is constructed, where x is an 8-dimensional homography parameter vector. The SVD (singular value decomposition) method is used to solve the optimal solution, and the scale consistency is maintained by matrix normalization. In order to improve the transformation accuracy, the matrix is pre-processed by coordinate normalization, and all corner point coordinates are recalibrated with the image center coordinate as the origin before solving, so that the transformation result is more stable. During the matrix calculation process, the reference coordinate system of the target image is fixed, and the rotation angle, scaling factor, translation vector and projection distortion component are extracted according to the spatial transformation behavior of the matching corner points, which are expressed as θ, s, t(x,y) and P respectively. By analyzing the asymmetric terms in the matrix, the degree of perspective projection distortion of the image is extracted. If the absolute value of the corresponding matrix elements H[2][0] and H[2][1] exceeds 0.002, it is considered that there is obvious projection distortion. All transformation parameters together with matrix elements are recorded in floating point form, and the complete transformation matrix data is output.

[0117] Step S313: calculating the stitching error of the image stitching distortion data according to the transformation matrix relationship;

[0118] In this embodiment, the coordinates of the corner points in the original image are transformed and mapped based on the homography transformation matrix identified in step S312 to obtain predicted transformed coordinates. The predicted coordinates are then compared with the observed coordinates of the corner points in the actual stitched image to calculate the stitching error for each point. The error is calculated using the Euclidean distance formula, with units of pixels. For all valid matching points, the error value and the corresponding image region position are recorded for each point, and the mean, variance, and maximum value of the error are calculated. The stitching error tolerance standard is set as: an average error of less than 1.2 pixels and a maximum error of no more than 2.5 pixels. If the average error of any image pair exceeds this standard, the stitching quality is marked as abnormal. In addition, the error points are projected back to the original image coordinate system to generate an error heat map. The error heat map uses bilinear interpolation to convert the error values into a color gradient image (red indicates high error areas, green indicates low error areas) and is output in pixel units to facilitate subsequent spatial analysis and geometric correction. Finally, the stitching error is output in a table format (including error value, corner point coordinates, and image number), accompanied by image visualization results.

[0119] Step S314: performing geometric correction based on the stitching error to obtain geometric distortion data.

[0120] In this embodiment, the stitching image is corrected using a local geometric deformation compensation method based on the stitching error calculated in step S313. The correction method is based on the Thin-Plate Spline (TPS) deformation model, which constructs a mapping relationship between the source point (original corner point coordinates) and the target point (ideal position coordinates) to generate a deformation function f(x, y). In actual operation, for each image with errors, the corner points with larger errors are used as control nodes (the error value exceeds 1.5 pixels), and the TPS model is input to remap the image coordinates. The image coordinate mapping uses bilinear interpolation to regenerate the corrected pixel grayscale value and reconstruct the entire image. Set the output image size to be the same as the original image size. Figure 1 To ensure the accuracy of the image quality, a boundary compensation strategy is implemented during the correction process, using a mirrored filling method to reduce interpolation errors. After correction, the image undergoes a new round of corner extraction and error verification to ensure that the geometric correction meets the error standard (the average error is controlled within 0.8 pixels). The corrected image data is ultimately output as the geometric distortion correction result. Together with the pre- and post-correction coordinate pairs and error comparison data, a geometric distortion dataset is formed, providing the foundation for geometric accuracy support for subsequent lens analysis and leaf area index estimation.

[0121] Preferably, step S32 is specifically as follows:

[0122] Step S321: calculating the image distortion based on the image stitching distortion data; and locating the lens shooting position according to the image distortion;

[0123] In this embodiment, in the image stitching distortion data, the image edge structure is extracted for each image separately, and the Canny edge detection algorithm is used, with the low threshold set to 50 and the high threshold set to 150. The detection results are subjected to Hough straight line transform to extract the main structural line segments in the image. The angular distribution and positional relationship of each line segment are recorded, and the root mean square error (RMSE) of all straight line segments in the image deviating from the ideal straight line is calculated. This value is set as the image distortion index, with the unit being pixel. Under the premise of knowing the fixed viewing angle and route arrangement structure of the camera carried by the drone, a corresponding relationship table between the image number and the GPS shooting coordinate is constructed. Combined with the spatial index of each image in the image stitching distortion data, images with a distortion value higher than 3 times the standard deviation of the average value of all images are selected and marked as abnormally distorted images. The latitude and longitude coordinates and flight altitude of the original shooting are located according to their numbers (obtained through the sensor height information in the flight control data), and the output is the lens shooting position data.

[0124] Step S322: performing scratch detection based on the lens shooting position to obtain lens scratch data;

[0125] In this embodiment, after the lens shooting position is determined, the optical system self-test information and the images recorded before and after the mission in the corresponding drone flight log are retrieved. The central area of the image (a square with the center of the image as the origin and a side length of 1 / 3 of the side length of the image) is selected as the main imaging area of the optical path, and the structural texture consistency detection is performed on this area. Using the local variance calculation method, the analysis window is set to 11×11 pixels, the sliding step is 3 pixels, the local standard deviation image is extracted, and abnormal high-frequency noise is detected. Linear texture distribution extraction is performed on the detected high-frequency area, and the directional gradient map is extracted by the Sobel operator. The area with a gradient direction aggregation distribution angle greater than 45 degrees is identified as the scratch area. The candidate scratch area is further subjected to image grayscale threshold segmentation, and the threshold is 60 units above the peak of the original image grayscale histogram. Morphological closing operation is used for processing, and the structural element adopts a 5×5 elliptical structure kernel to obtain a complete scratch contour area. The total length, maximum length, average width, and density (scratch length per square centimeter) of the plan are calculated and bound to the image number and shooting time to form a lens scratch data table. All values are converted to standard millimeters. Lengths are calculated based on image size and focal length calibration parameters, which are provided in the camera's factory documentation.

[0126] Step S323: Acquire wind data and perform data preprocessing to obtain wind data to be analyzed;

[0127] In this embodiment, during the stage of acquiring wind data, the data from the on-site unmanned meteorological observation module is accessed. The data format is timestamp + wind speed (m / s) + wind direction (0-360 degrees) + mutation frequency, the sampling frequency is 1Hz, and the data source is data continuously recorded for 10 minutes during the flight. First, data denoising is performed, and a sliding mean filter (window length is 5 seconds) is used to smooth the wind speed and wind direction data sequences. For the sudden wind speed point, the wind speed change is greater than 2.5m / s within 5 seconds as the wind speed mutation point, and the interpolation method (linear interpolation) is used to complete it. The wind speed vector is projected to the flight direction of the drone to obtain the relative wind speed component. The flight direction is inversely calculated based on the coordinate difference between two consecutive GPS points, and the angle between the wind speed and the flight direction is calculated. The impact wind speed component along the flight direction is extracted by the cosine theorem. Finally, the wind speed data is matched according to the shooting time of each image to construct the image number-wind speed-wind direction triple data to form the wind data to be analyzed. The storage format is CSV format, the time accuracy is unified to millisecond level, and the unit is strictly maintained in m / s and angle (°).

[0128] Step S324: performing wind impact simulation on the lens scratch data according to the wind data to be analyzed to obtain wind impact data;

[0129] In this embodiment, the wind data to be analyzed obtained in step S323 is used as input parameters to model the camera lens structure corresponding to the lens scratch data in step S322. The outer surface diameter of the lens is defined as 25 mm, the radius of curvature is 12.5 mm, the material is high-transmittance glass, and the density is 2.6 g / cm 3 The surface roughness is set to Ra = 0.001μm. The wind impact simulation uses a CFD two-dimensional simulation method based on the finite difference method to load the lens surface with wind pressure per unit area. At each wind time point, the direction component of the force is generated based on the wind speed and wind direction. The wind pressure is calculated according to the standard formula P = 0.5×ρ×v 2 Calculation is performed, where the air density ρ is taken as 1.225 kg / m 3 , v is the component velocity along the normal direction of the lens. The wind pressure is applied to the identified scratch area, and local wind pressure disturbances are applied to different areas according to the spatial distribution of the scratches. The simulation area is divided into 0.2mm×0.2mm grids, and an explicit time step of 0.001 seconds is used for 10 seconds. The force history sequence on each scratch area is obtained by simulation, and the total pressure per unit area is output in N / m 2 , and mark the pressure increase point caused by the sudden change of wind speed (defined as the force difference between two adjacent moments is greater than 5N / m 2 ). Finally, a wind impact data file is formed, which includes the time point, scratch number, force value and wind speed direction angle.

[0130] Step S325: identifying the movement trajectory of sand particles based on the wind impact data;

[0131] In this embodiment, based on the wind impact data, a model of the initial state of sand particles in the wind field is constructed. Assuming that the density of sand particles is 2650 kg / m 3 The particle size is 0.3 mm and the initial height is 30 m. The Lagrangian method is used to simulate the trajectory of sand particles. The initial velocity is the direction vector of the wind speed impact value. The time step is set to 0.01 seconds. The total simulation time is 10 seconds. The gravity acceleration is 9.81 m / s. 2 The simulation of sand particle trajectory takes into account the resistance model F_d=0.5×C_d×A×ρ×v 2 , the drag coefficient C_d is taken as 0.47, and the cross-sectional area A is calculated as πd based on the spherical body 2 / 4. Wind direction data is dynamically loaded as a two-dimensional vector, and wind speed is updated once per second. During the simulation, the spatial coordinates (x, y, z) of each sand particle's trajectory are recorded to determine whether it impacts the drone's camera during flight. The judgment criteria are: the distance between the sand particle and the camera's center is less than 1 cm, and the error between the sand particle's height and the camera's horizontal level is less than 1.5 cm. The output trajectory data consists of the sand particle number, time point, three-dimensional coordinates, and impact status (0 / 1). The trajectories of all sand particles that successfully impact are counted and their numbers, trajectories, and impact angles are recorded as the final sand particle trajectory data.

[0132] Step S326: performing impact force detection on the lens scratch data according to the sand movement trajectory to obtain the sand impact force;

[0133] In this embodiment, the impact point in the sand trajectory output in step S325 is used to extract the end value of the sand velocity vector and calculate the kinetic energy E_k=0.5×m×v based on the mass of the sand. 2 , in joules. Then, based on the impact angle θ (the angle with the normal of the lens surface), calculate the component F in the actual force direction = E_k × cos(θ) / Δs, where Δs is set to the action time 0.001 seconds multiplied by the impact speed, to obtain the impact force per unit time, in Newtons. Perform coordinate matching on the scratch area where the impact occurs, map the impact force value to the existing scratch area number, and update the cumulative number of impacts and the maximum single impact force of each scratch area. Count the total impact force of each lens area (summed by time integral) in N·s to form a sand impact force data table. The data is composed of image number, scratch number, total number of impacts, maximum impact value and total impact force, accurate to three decimal places.

[0134] Step S327: predicting lens damage at the lens shooting position based on the impact force of the sand particles to obtain lens damage data.

[0135] In this embodiment, combined with the lens material parameters (glass impact strength is set to 15N / mm 2 ), calculate the impact force density for each scratch area, normalize it by area, and convert the unit to N / mm 2 The judgment standard is: the impact density exceeds 5N / mm 2 The scratched area enters a high-risk damage state when the pressure exceeds 10N / mm 2When the lens is damaged, it is considered to be in a critical damage state. All scratched areas are grouped according to image number, and images with high-risk and critical areas are marked as damage warning images. For all damaged images, the lens shooting location (latitude and longitude, flight altitude), time, and damage level (0-damageless, 1-high risk, 2-critical) are summarized to form a lens damage data table with data accuracy to five decimal places. The final output is a structured data file with fields including image number, location coordinates, damage level, maximum impact force value, maximum impact density, scratch number and timestamp.

[0136] Preferably, step S4 is specifically as follows:

[0137] Step S41: capturing a high-precision orthophoto of the terrain-undulating area according to the adaptively adjusted focal length data;

[0138] In this embodiment, a high-resolution camera mounted on a drone is used to capture images of rugged terrain. To ensure geometric accuracy of the images, the focal length is adaptively optimized using an automatic adjustment algorithm. The focal length is set between 18mm and 55mm. The camera focal length is adjusted based on the complexity of the terrain, the shooting altitude, and the required ground pixel resolution to ensure sufficient detail and clarity in the captured images. During flight, the flight control system monitors the flight altitude, attitude angle, and GPS positioning information in real time. Combining the flight trajectory with ground feature data in the shooting area, the camera focal length is dynamically adjusted to ensure optimal viewing angle and clarity for each frame. During the shooting process, the overlap of each image is set to 60% to ensure sufficient coverage and provide reference feature points for subsequent stitching. Pre-flight exposure parameters, white balance, and ISO sensitivity for all images are uniformly set to 1 / 1000 second, 200 ISO, and daylight mode for white balance. During the shooting process, data from each automatic camera focus adjustment is collected, and the focal length value and corresponding flight altitude are recorded for subsequent image stitching and orthophoto generation. Using data and images transmitted by drones, adaptive image processing technology generates high-precision orthophotos of terrain-rough areas. These orthophotos are integrated using image stitching algorithms to minimize distortion of the imaged area and are calibrated to the geographic coordinate system.

[0139] Step S42: Repairing the distortion of the orthophoto of the terrain undulating area according to the high-precision orthophoto to obtain a repaired orthophoto;

[0140] In this embodiment, the generated high-precision orthophoto is subjected to distortion analysis. Image registration technology is used to align images taken from different angles. The core parameters of image registration include extracting key points in the image using a feature matching algorithm (such as SIFT or ORB algorithm). These key points need to correspond to the same feature position in multiple images. By calculating the transformation matrix between images, geometric distortion caused by factors such as flight trajectory, perspective difference, and camera lens distortion is eliminated. When repairing the distorted image, the perspective correction of different image blocks is performed using the projection transformation method, and the correction transformation matrix of each image is calculated using the least squares method to ensure accurate docking of the image coordinate system. In addition, the distortion tolerance of the repair algorithm is set to 0.5 pixels, that is, when the deformation error after image repair is less than 0.5 pixels, the part of the image is considered to have been successfully repaired. Bilinear interpolation is used to fill image pixels to ensure that the visual quality of the image is maintained after repair. After the image repair is completed, the repaired orthophoto is spatially calibrated using a geographic information system (GIS) to ensure that each pixel has an accurate geographic coordinate position. During this process, geographic coordinate support is provided by GPS and IMU data to complete the geographic registration of the orthophoto, and finally a high-precision orthophoto is obtained after distortion repair.

[0141] It is particularly important that step S42 includes the following steps:

[0142] Step S421: performing spectral correction on the orthophoto of the terrain undulating area according to the high-precision orthophoto to obtain spectral correction data;

[0143] In this embodiment, high-precision orthophoto data is acquired and spectrally corrected. The purpose of spectral correction is to correct for deviations in the image's spectral values caused by factors such as atmosphere, moisture, and illumination. This process utilizes radiometric correction techniques, primarily consisting of atmospheric correction and spectral restoration. First, atmospheric correction is performed on the image based on a known standard atmospheric model (such as the MODTRAN model). This correction takes into account the absorption and scattering effects of different bands in the atmosphere and corrects the reflectance values of each band. The specific steps of atmospheric correction include inputting the image's sensor data, selecting an appropriate radiation transfer model, and calculating and applying correction factors based on parameters such as the image's capture time, geographic location, and climatic conditions. Then, spectrally restoration is performed on the corrected image based on known ground features or reflectance standards (such as reflectance standards for land, vegetation, and water). During the restoration process, spectral comparison and adjustment are performed using ground calibration points. Through linear or nonlinear transformations of multi-band data, inconsistent spectral responses are compensated to ensure that the reflectance values of all bands are within the standard range. Key parameters in this process include band response correction coefficients, ground feature reflectance standards, and meteorological data. Through these steps, the corrected spectral data is finally obtained, providing accurate spectral information for subsequent image analysis.

[0144] Step S422: performing geometric distortion repair on the orthophoto of the terrain undulating area according to the high-precision orthophoto to obtain geometric distortion repair data;

[0145] In this embodiment, high-precision orthophotos are used to analyze images of areas with undulating terrain to identify areas of geometric distortion in the image. Geometric distortion is typically caused by deviations in the image's geometric shape due to factors such as shooting angle, terrain undulation, and sensor error. First, image registration techniques are used to identify distorted areas using control points with known geographic coordinates. Image registration determines the location and magnitude of image distortion by matching control points in the image with known standard geographic coordinates. Next, image transformation techniques, such as backprojection and least-squares optimization, are used to geometrically correct the distorted areas in the image. During this process, the image coordinate system is converted to the geographic coordinate system to restore the geometry of the distorted areas. In specific implementation, a high-precision reference coordinate system, such as WGS84 or UTM, is selected to align the pixel coordinates of the distorted image with the actual geographic coordinates. During the restoration process, bilinear interpolation is used to resample the pixels in the restoration area, and the transformation parameters are optimized and adjusted based on the control point errors. For areas with severe distortion, a nonlinear optimization algorithm is used for depth correction. Finally, through this process, the image data after geometric distortion repair is obtained, which provides an accurate geometric reference for further data analysis.

[0146] Step S423: Integrate the spectral correction data and the geometric distortion restoration data to obtain a restored orthophoto.

[0147] In this embodiment, after obtaining the spectral correction data and the geometric distortion repair data, the last step is to integrate the two to obtain the repaired orthophoto. The integration process first involves data fusion of the completed spectral correction data and geometric repair data to ensure that the two data sources have the same spatial resolution and geographic reference. The specific steps include aligning the two data through an image registration algorithm to ensure that their spatial positions are completely consistent. At this time, if there are differences in the coordinate systems of the two, the coordinate systems need to be unified, and a high-precision geographic coordinate system (such as WGS84) is usually used for conversion to ensure that all image data have a consistent coordinate system and spatial resolution. After completing the coordinate alignment, the spectrally corrected image and the geometrically repaired image are merged into a whole image through pixel-level fusion technology (such as multi-band splicing or resampling technology). This process requires the selection of an appropriate resampling method (such as nearest neighbor interpolation or bilinear interpolation) to avoid information loss or distortion during the image fusion process. At the same time, when fusing images, it is necessary to ensure the accuracy of the data to avoid inconsistencies between the two data affecting the quality of the final image. Finally, the generated restored orthophoto data will serve as the basic data for subsequent analysis (such as crop leaf area index (LAI) estimation), providing reliable image data support for higher-precision crop monitoring and analysis.

[0148] Step S43: performing plant species identification based on the restored orthophoto to obtain plant species data;

[0149] In this embodiment, plant species identification is performed using a deep learning image classification method using restored high-precision orthophotos. First, a convolutional neural network (CNN) algorithm is used to automatically identify plant areas in orthophotos. The dataset used for training contains high-resolution images of different plant species, with the species of each plant labeled. The image preprocessing steps include cropping and scaling. All input images are adjusted to 256×256 pixels, and data enhancement (such as rotation, translation, and brightness changes) is performed to improve the robustness of recognition. Plant species classification is performed using a pre-trained deep learning model. The model outputs a probability distribution of plant species categories, and the threshold is set to 0.8, that is, when the recognition probability of a species is greater than 80%, the area is considered to be that species. In order to further optimize the recognition results, a clustering algorithm (such as K-means) is used to spatially aggregate similar species to reduce the probability of misidentification. The plant species data includes the species name, location coordinates, and confidence level of each identified area. Plant species data are recorded in the database and spatially located in combination with the coordinate information of the image. The data table contains fields such as species category, area, longitude and latitude, and identification confidence.

[0150] Step S44: estimating the crop leaf area index based on the plant species data and the plant growth stage.

[0151] In this embodiment, the growth stage of each crop is further obtained in combination with the identified plant species data. The timestamp at the time of flight is compared with the growth pattern of the crop species, and the current growth stage of the crop is determined in combination with meteorological data (temperature, precipitation, light intensity, etc.). For each crop, a standard time node of the growth stage is set, for example: 30 days after sowing is the seedling stage, 60 days is the growth period, and 120 days is the maturity period. The crop leaf area index (LAI) is estimated by the relationship formula between the vegetation index (NDVI, Normalized Difference Vegetation Index) and the crop growth stage. NDVI is calculated by the following formula:

[0152]

[0153] NIR represents the reflectance of the near-infrared band, and Red represents the reflectance of the red band. NDVI values range from -1 to 1, with higher values indicating denser vegetation cover.

[0154] The leaf area index (LAI) is calculated using a polynomial regression model based on the empirical regression relationship between NDVI and LAI. Different regression formulas are set for different crop growth stages. For example, a linear relationship between LAI and NDVI is used in the early growth stage, and a logarithmic relationship between LAI and NDVI is used in the mature stage. The formula is:

[0155] LAI = a × (NDVI) b ;

[0156] Where a and b are regression coefficients obtained by fitting historical data based on crop type and growth stage. The leaf area index value of each plant species is finally obtained and recorded in the crop leaf area index data table by crop type, location, and leaf area index. The data unit is square meters per square meter.

[0157] The present invention is therefore intended to be illustrative and non-restrictive in all respects, with the scope of the invention being defined by the appended claims rather than the foregoing description, and all changes that come within the meaning and range of equivalents of the application documents are intended to be embraced therein.

[0158] The foregoing description is intended only to provide specific embodiments of the present invention, which will enable those skilled in the art to understand and implement the present invention. Various modifications to these embodiments will be readily apparent to those skilled in the art, and the general principles defined herein may be implemented in other embodiments without departing from the spirit or scope of the present invention. Therefore, the present invention is not intended to be limited to the embodiments shown herein, but is to be construed in the widest possible manner consistent with the principles and novel features disclosed herein.

Claims

1. A crop leaf area index estimation method based on UAV image processing, characterized in that: The following steps are involved: Step S1: obtaining a remote sensing image of the area monitored by the drone; identifying the terrain undulating area based on the remote sensing image of the area monitored by the drone, and reconstructing the orthophoto of the terrain undulating area to generate an orthophoto of the terrain undulating area; Step S2: identifying terrain-blocked areas based on remote sensing images of the area monitored by the drone; performing image stitching distortion detection on the terrain-blocked areas to obtain image stitching distortion data; performing spectral interpolation on the gap areas based on the image stitching distortion data, and determining the plant growth stage in the gap areas after spectral interpolation; Step S3: performing lens distortion detection based on the image stitching distortion data to obtain lens distortion data; Identify focus component damage characteristics based on lens distortion data; Adaptively adjust the focal length based on the damage characteristics of the focusing component to obtain adaptive adjustment focal length data; Step S4: performing plant species identification on the orthophoto of the terrain undulating area according to the adaptively adjusted focal length data to obtain plant species data; and estimating the crop leaf area index according to the plant species data and the plant growth stage.

2. The crop leaf area index estimation method based on drone image processing according to claim 1 is characterized in that: Step S1 is specifically as follows: Step S11: Acquire remote sensing images of the drone monitoring area; Step S12: generating point cloud data of the monitoring area according to the remote sensing image of the monitoring area by the UAV, and generating a digital elevation model based on the point cloud data of the monitoring area; Step S13: extracting regional pixel height values according to the digital elevation model, and calculating the slope of the monitoring area based on the regional pixel height values; Estimate spectral reflectance values from a digital elevation model and identify shadow areas based on the spectral reflectance values; Step S14: determining the terrain undulation area based on the slope of the monitoring area and the shadow area; Step S15: reconstructing the orthophoto of the terrain undulating area to generate an orthophoto of the terrain undulating area.

3. The crop leaf area index estimation method based on drone image processing according to claim 2 is characterized in that: Step S15 is specifically as follows: Step S151: collecting remote sensing images of the terrain undulating area to obtain remote sensing images of the terrain undulating area; Step S152: Acquire ground control point data; extract corner points of the remote sensing image of the terrain relief area, wherein the number of extracted corner points is set to 50-100; perform local feature point registration based on the ground control point data and the corner points to obtain local registration data; Step S153: performing image projection conversion on the orthophoto of the terrain relief area according to the local registration data to obtain a local orthophoto projection; Step S154: performing texture correction on the local orthographic projection to obtain a texture-corrected local orthographic projection; Step S155: performing image stitching based on the texture-corrected local orthophoto projection to generate an orthophoto of the terrain relief area.

4. The crop leaf area index estimation method based on UAV image processing according to claim 1, characterized in that: The identification of terrain-blocked areas in step S2 is specifically as follows: Grayscale conversion is performed based on the remote sensing image of the UAV monitoring area to obtain a grayscale remote sensing image; Calculate the gray-level co-occurrence matrix of gray-level remote sensing images; The roughness is calculated based on the gray-level co-occurrence matrix; According to the gray-level co-occurrence matrix, the entropy value is calculated; Identify uneven texture areas in grayscale remote sensing images based on roughness and entropy values; Detect linear transition edges based on remote sensing images of the area monitored by UAVs; Identify ridge areas in uneven texture areas based on linear transition edges; obtain sunlight angles; perform lighting simulation on ridge areas based on sunlight angles, and identify terrain occlusion areas during the simulation process.

5. The crop leaf area index estimation method based on UAV image processing according to claim 1, characterized in that: The identification of terrain-blocked areas in step S2 is specifically as follows: Collect multiple remote sensing images of terrain-blocked areas to obtain multiple remote sensing images; Generate terrain occlusion stitching images based on multiple remote sensing images; identify terrain occlusion edge segments based on terrain occlusion stitching images; Calculating the line segment slope of the terrain occlusion edge segment; identifying the angle jump line segment of the terrain occlusion edge segment based on the line segment slope; Extract terrain occlusion connected areas based on terrain occlusion stitching images; Find 8 adjacent edge pixels based on the terrain occluded connected area, and expand the edge chain based on the 8 adjacent edge pixels to obtain edge chain data; Detect edge chain breakpoints based on edge chain data; Determine the geometric structure distortion area image of the terrain occlusion stitching image based on the angle jump line segment and the edge chain break point; Extract image acquisition timestamp based on terrain occlusion stitching image; obtain UAV three-axis angle data; Map the image acquisition timestamp to the drone's three-axis angle data to obtain the image acquisition drone's three-axis angle data; Based on the image acquisition of the UAV's three-axis angle data, the UAV attitude jump detection is performed to obtain the UAV attitude jump data; Extracting attitude jump image frames based on the attitude jump data of the UAV; mapping the attitude jump image frames to the terrain occlusion stitching image to obtain the attitude jump terrain occlusion image; The posture jump terrain occlusion image and the geometric structure distortion area image are integrated to obtain the image stitching distortion data.

6. The crop leaf area index estimation method based on drone image processing according to claim 1, characterized in that: The spectral interpolation of the gap region in step S2 is specifically as follows: Extract zero pixel data based on image stitching distortion data; locate gap areas based on zero pixel data; A 5x5 neighborhood window is used to extract the effective spectrum value of the gap area; Perform weighted average calculation based on the effective spectrum value to obtain a weighted spectrum value; Performing spectral interpolation on the gap region based on the weighted spectral value to obtain a spectral interpolation gap region; Calculating vegetation index based on spectral interpolation gap area; identifying vegetation area in spectral interpolation gap area according to vegetation index; Collect leaf spectra of vegetation areas and invert chlorophyll content based on the leaf spectra; Determine plant growth stage based on chlorophyll content.

7. The crop leaf area index estimation method based on UAV image processing according to claim 1, characterized in that: Step S3 is specifically as follows: Step S31: performing geometric distortion identification based on the image stitching distortion data to obtain geometric distortion data; Step S32: performing lens damage detection based on the image stitching distortion data to obtain lens damage data; Step S33: Integrate the geometric distortion data and the lens damage data to obtain lens distortion data; Step S34: identifying focus assembly damage features based on the lens distortion data; Step S35: performing adaptive focal length adjustment based on the damage characteristics of the focusing component to obtain adaptively adjusted focal length data.

8. The crop leaf area index estimation method based on UAV image processing according to claim 7, characterized in that: Step S31 is specifically as follows: Step S311: extracting image corner points based on the image stitching distortion data, and performing corner point matching based on the image corner points to obtain matching corner point data; Step S312: identifying a transformation matrix relationship based on the matching corner point data; Step S313: calculating the stitching error of the image stitching distortion data according to the transformation matrix relationship; Step S314: performing geometric correction based on the stitching error to obtain geometric distortion data.

9. The crop leaf area index estimation method based on UAV image processing according to claim 7, characterized in that: Step S32 is specifically as follows: Step S321: calculating the image distortion based on the image stitching distortion data; and locating the lens shooting position according to the image distortion; Step S322: performing scratch detection based on the lens shooting position to obtain lens scratch data; Step S323: Acquire wind data and perform data preprocessing to obtain wind data to be analyzed; Step S324: performing wind impact simulation on the lens scratch data according to the wind data to be analyzed to obtain wind impact data; Step S325: identifying the movement trajectory of sand particles based on the wind impact data; Step S326: performing impact force detection on the lens scratch data according to the sand movement trajectory to obtain the sand impact force; Step S327: predicting lens damage at the lens shooting position based on the impact force of the sand particles to obtain lens damage data.

10. The crop leaf area index estimation method based on UAV image processing according to claim 1, characterized in that: Step S4 is specifically as follows: Step S41: capturing a high-precision orthophoto of the terrain-undulating area according to the adaptively adjusted focal length data; Step S42: Repairing the distortion of the orthophoto of the terrain undulating area according to the high-precision orthophoto to obtain a repaired orthophoto; Step S43: performing plant species identification based on the restored orthophoto to obtain plant species data; Step S44: estimating the crop leaf area index based on the plant species data and the plant growth stage.

Images