A method for radiometric consistency correction of aircraft multispectral images

Abstract

The invention relates to a method for correcting radiation consistency of aircraft multispectral images. The method comprises the following steps: firstly, performing vignetting correction; then performing histogram contrast stretching on the image; then performing feature matching on a reference image and an image after histogram contrast stretching by using a SIFT operator to obtain matching point pairs; then performing linear regression based on the pixel values ​​of the matching point pairs according to a RANSAC algorithm to establish a linear correction model; finally, applying the linear correction model, performing radiation consistency correction on the image to be corrected, and performing bilateral filtering on the corrected image to remove noise points; the design scheme overcomes the drawback that the traditional correction method cannot effectively obtain sufficient same-name points in low-illuminance images, and thus cannot carry out radiation information correction between images; at the same time, the modeling method has strong robustness, effectively resists interference from noise points in the image, can make the radiation conditions of thousands of unmanned aerial vehicle images in one sortie tend to be consistent, and eliminates the image distortion problem caused by strong reflection of ground objects and drastic changes in light intensity; thereby providing solid data support for subsequent quantitative remote sensing research.

Description

Technical Field

[0001] The invention relates to a method for correcting the radiation consistency of aircraft multi-spectral images, and belongs to the technical field of unmanned aerial vehicle remote sensing. Background Art

[0002] Multispectral drones play an important role in the study of crop moisture, growth, nitrogen content, and pest and disease monitoring. However, the data obtained by drones often encounter some problems due to environmental factors, such as uneven brightness of the same flight image due to drastic changes in light intensity in cloudy weather. Even in clear and cloudless weather, when there are reflective objects such as water, glass, and stainless steel on the ground, due to different shooting angles, some images will produce flares, and the overall darkening of areas other than the reflective areas will appear. This will cause the same object to present serious spectral feature inconsistencies in different images, and the spectral distortion of the object will appear in the subsequent drone image stitching results.

[0003] In order to improve the accuracy of subsequent quantitative research, it is necessary to perform radiation consistency correction on the images collected by the multispectral camera of the UAV. At present, the radiation consistency correction between remote sensing images mainly adopts the histogram matching method, the correction method based on real-time illumination information, and the correction method based on statistical regression. The above methods often achieve relatively good results in correcting the inconsistent radiation information between images of the same flight caused by changes in illumination. However, when there are low-illuminance images caused by reflections from ground objects, due to the low overall contrast of the image, few visual details, and low average brightness, it is impossible to select enough high-quality sample sets manually or by algorithm, so the above methods are usually not ideal when processing such images. Summary of the invention

[0004] The technical problem to be solved by the present invention is to provide a method for correcting the radiation consistency of aircraft multispectral images. The method adopts a new design method, can overcome variable scenes, obtain reliable correction results, and effectively correct the radiation information between images while removing noise points in the image as much as possible.

[0005] In order to solve the above technical problems, the present invention adopts the following technical solutions: the present invention designs a method for correcting the radiation consistency of aircraft multispectral images, based on the multispectral reference images obtained by the aircraft shooting the target type area at each preset flight altitude, for each multispectral image to be corrected of the target type area shot by the aircraft at the same flight as the multispectral reference image and at each preset flight altitude, the radiation consistency correction is realized; the correction method comprises executing the following steps A to G for each multispectral image to be corrected;

[0006] Step A. Apply the radial vignetting correction model to perform vignetting correction on the multispectral image to be corrected, update the multispectral image to be corrected, and then proceed to step B;

[0007] Step B. Performing histogram contrast stretching processing on the multispectral image to be corrected to obtain the multispectral stretched image to be corrected, and then proceeding to step C;

[0008] Step C. identifying each feature point in the multispectral stretched image to be corrected, and extracting the feature vectors corresponding to each feature point, and then proceeding to step D;

[0009] Step D. Based on the feature points in the multispectral reference image with the same flight altitude as the multispectral image to be corrected obtained in advance according to the method in step C, and the feature vectors corresponding to the feature points, by comparing the Euclidean distances between the feature points, obtain the feature point groups that match each other in the multispectral stretched image to be corrected and the multispectral reference image, and then obtain the corresponding feature point position groups, and then enter step E;

[0010] Step E. For each feature point position group, respectively, obtain the pixel value of the corresponding position in the feature point position group in the multispectral reference image with the same flight altitude as the multispectral image to be corrected, as the reference pixel value, and obtain the pixel value of the corresponding position in the feature point position group in the multispectral image to be corrected, as the pixel value to be corrected, and form a pixel fitting sample group corresponding to the feature point position group by combining the pixel value to be corrected and the reference pixel value; then obtain the pixel fitting sample groups corresponding to each feature point position group, and then enter step F;

[0011] Step F. Based on the correspondence between the pixel values ​​to be corrected and the reference pixel values ​​in each pixel fitting sample group, data fitting is performed for a preset model according to a preset fitting method to form an image radiation consistency correction model corresponding to the flight altitude of the aircraft and the multispectral image to be corrected, and then proceeding to step G;

[0012] Step G. Apply the image radiation consistency correction model to correct and update the pixel values ​​at each position in the multispectral image to be corrected, obtain the multispectral correction image corresponding to the multispectral image to be corrected, and then enter step H;

[0013] Step H: performing bilateral filtering on the multi-spectral correction image to update the multi-spectral correction image.

[0014] As a preferred technical solution of the present invention: a calibration plate is placed on the ground in an area of ​​the target type area that is flat, far away from water and has no shadows, and the aircraft performs multispectral photography of the area at each preset flight altitude to obtain multispectral reference images of the target type area corresponding to each preset flight altitude, and for each multispectral reference image, a radial vignetting correction model is applied to perform vignetting correction on the multispectral reference image, and each multispectral reference image is updated.

[0015] As a preferred technical solution of the present invention: the step B includes steps B1 to B3 as follows:

[0016] Step B1. Count the number of pixels corresponding to each different pixel value in the multispectral image to be corrected, and form the pixel ratio corresponding to each different pixel value through the ratio of the number to the total number of pixels in the multispectral image to be corrected, and then obtain the sum of the pixel ratios corresponding to each pixel value not greater than the pixel value for each different pixel value as the cumulative pixel ratio corresponding to the pixel value, and then obtain the cumulative pixel ratio corresponding to each different pixel value, and then construct a histogram with the horizontal axis as the pixel value and the vertical axis as the cumulative pixel ratio, and enter step B2;

[0017] Step B2. Based on the cumulative pixel percentages in the histogram, a cutoff value comparison method is applied to obtain the lower pixel value and the upper pixel value of each different pixel value in the multispectral image to be corrected, and then proceed to step B3;

[0018] Step B3. Use the lower pixel value and the upper pixel value to replace the maximum pixel value and the minimum pixel value in the formula involved in contrast stretching of the multispectral image to be corrected, and then perform contrast stretching on the multispectral image to be corrected to obtain the multispectral stretched image to be corrected, and then enter step C.

[0019] As a preferred technical solution of the present invention: the step B2 is the following step B2-I;

[0020] Step B2-I. Based on the cumulative pixel ratios in the histogram, combined with the preset ratio cutoff value a and b=100%-a, determine whether there is a cumulative pixel ratio equal to a. If so, use the cumulative pixel ratio as the target cumulative pixel ratio, otherwise select the cumulative pixel ratio closest to a as the target cumulative pixel ratio; at the same time, determine whether there is a cumulative pixel ratio equal to b. If so, use the cumulative pixel ratio as the target cumulative pixel ratio, otherwise select the cumulative pixel ratio closest to b as the target cumulative pixel ratio; then obtain the pixel values ​​corresponding to the two target cumulative pixel ratios, use the smaller pixel value of the two pixel values ​​as the lower pixel value and the larger pixel value as the upper pixel value, and then enter step B3.

[0021] As a preferred technical solution of the present invention: in step B2-I, when the preset percentage cutoff value a is less than 50%:

[0022] If the number of cumulative pixel ratios closest to a is 2, the smaller cumulative pixel ratio of the two cumulative pixel ratios is selected as the target cumulative pixel ratio, and the pixel value corresponding to the target cumulative pixel ratio is obtained as the next pixel value;

[0023] If the number of cumulative pixel ratios closest to b is 2, the larger cumulative pixel ratio of the two cumulative pixel ratios is selected as the target cumulative pixel ratio, and the pixel value corresponding to the target cumulative pixel ratio is obtained as the upper pixel value;

[0024] When the preset percentage cutoff value a is greater than 50%:

[0025] If the number of cumulative pixel ratios closest to a is 2, the larger cumulative pixel ratio of the two cumulative pixel ratios is selected as the target cumulative pixel ratio, and the pixel value corresponding to the target cumulative pixel ratio is obtained as the upper pixel value;

[0026] If the number of cumulative pixel point ratios closest to b is 2, the smaller cumulative pixel point ratio of the two cumulative pixel point ratios is selected as the target cumulative pixel point ratio, and the pixel value corresponding to the target cumulative pixel point ratio is obtained as the next pixel value.

[0027] As a preferred technical solution of the present invention: the step B2 is the following step B2-II;

[0028] Step B2-II. Based on the cumulative pixel percentages in the histogram sorted from small to large, the number N of cumulative pixel percentages is obtained, combined with the preset percentage threshold c of less than 50%, according to Select the dth cumulative pixel point ratio in sequence as the target cumulative pixel point ratio, and select the dth cumulative pixel point ratio from the bottom as the target cumulative pixel point ratio, and obtain the pixel values ​​corresponding to the two target cumulative pixel point ratios respectively, use the smaller pixel value of the two pixel values ​​as the lower pixel value and the larger pixel value as the upper pixel value, and then enter step B3.

[0029] As a preferred technical solution of the present invention: in the step B3, the lower pixel value and the upper pixel value are used to replace the maximum pixel value b and the minimum pixel value a in the following formula for contrast stretching of the multispectral image to be corrected, respectively.

[0030]

[0031] Then, the pixel value of each pixel point in the multispectral image to be corrected is updated by the above formula to realize contrast stretching of the multispectral image to be corrected, and obtain the multispectral stretched image to be corrected, wherein f(x, y) represents the pixel value of the pixel point at position (x, y) in the multispectral image to be corrected, g(x, y) represents the pixel value after the pixel value of the pixel point at position (x, y) in the multispectral image to be corrected is updated, and c and d are 0 and 255 respectively.

[0032] As a preferred technical solution of the present invention: SIFT algorithm is applied in step C, and steps C1 to C5 are performed as follows:

[0033] Step C1. Based on the pixel values ​​of the pixels in the image, a Gaussian function is applied to blur and downsample the multispectral stretched image to be corrected, and an image Gaussian pyramid is constructed, and then step C2 is entered;

[0034] Step C2. In the obtained Gaussian difference pyramid space, for each pixel point, the pixel point is compared with its 8 adjacent pixel points and 2×9=18 points in the adjacent upper and lower layers, a total of 26 points. If the pixel point is the maximum or minimum, it is used as a local key point; then each local key point is obtained, and then step C3 is entered;

[0035] Step C3. Based on each local key point, the position and scale of each local key point are accurately determined by fitting a three-dimensional quadratic function, and the edge response point is detected according to the Harris algorithm, and removed from all local key points, thereby obtaining the feature points in the multispectral stretched image to be corrected, and then entering step C4;

[0036] Step C4. For each feature point in the multispectral stretched image to be corrected, obtain the gradient between the feature point and its adjacent pixel points in each direction, and select the direction of the pixel point corresponding to the maximum gradient as the direction corresponding to the feature point; then obtain the directions corresponding to each feature point, and then enter step C5;

[0037] Step C5. For each feature point in the multispectral stretched image to be corrected, a feature vector corresponding to the feature point is constructed using the direction corresponding to the feature point, the position of the feature point in the image, and the gradient between the feature point and its adjacent pixel points in each direction; and then the feature vectors corresponding to each feature point in the multispectral stretched image to be corrected are obtained.

[0038] As a preferred technical solution of the present invention: in the step D, based on each feature point in the multispectral reference image with the same flight altitude as the multispectral image to be corrected obtained in advance according to the method of step C, and the feature vectors corresponding to each feature point, for each feature point in the multispectral stretched image to be corrected, the feature point is used as a candidate feature point, and the following steps D1 to D2 are performed to search for feature points matching the candidate feature points in the multispectral reference image, obtain each feature point group matching each other in the multispectral stretched image to be corrected and the multispectral reference image, and each corresponding feature point position group, and then enter step E;

[0039] Step D1. Obtain the Euclidean distance between each feature point in the multispectral reference image by obtaining the Euclidean distance between each feature point based on the feature vector corresponding to the feature point, and obtain the ratio between the closest distance and the second closest distance, and then proceed to step D2;

[0040] Step D2. Determine whether the ratio is less than a preset distance ratio threshold. If so, determine that the feature point to be selected matches the feature point in the multispectral reference image corresponding to the nearest distance to each other to form a group of feature points, and obtain the corresponding position of the group of feature points; otherwise, determine that there is no feature point in the multispectral reference image that matches the feature point to be selected.

[0041] As a preferred technical solution of the present invention: in the step F, based on the correspondence between the pixel values ​​to be corrected and the reference pixel values ​​in each pixel fitting sample group, data fitting is performed for a preset model according to the random sampling consistency algorithm RANSAC to form an image radiation consistency correction model for the flight altitude of the aircraft corresponding to the multispectral image to be corrected.

[0042] The method for correcting the radiation consistency of aircraft multispectral images described in the present invention has the following technical effects compared with the prior art by using the above technical solution:

[0043] (1) The method for correcting the radiation consistency of aircraft multispectral images designed by the present invention first performs vignetting correction; then performs histogram contrast stretching on the image; then uses the SIFT (Scale invariant feature transform) operator to perform feature matching on the reference image and the image after histogram contrast stretching to obtain matching point pairs; then based on the pixel values ​​of the matching point pairs, linear regression is performed according to the RANSAC (Random sample consensus) algorithm to establish a linear correction model; finally, the linear correction model is applied to perform radiation consistency correction on the image to be corrected, and bilateral filtering is performed on the corrected image to remove noise points; the design scheme overcomes the disadvantage that the traditional correction method cannot effectively obtain sufficient same-name points in low-illuminance images and cannot carry out radiation information correction between images; at the same time, the modeling method has strong robustness and effectively resists the interference of noise points in the image, and can make the radiation conditions of thousands of drone images in a sortie tend to be consistent, eliminating the image distortion problem caused by strong reflection of ground objects and drastic changes in light intensity; thereby providing solid data support for subsequent remote sensing quantitative research;

[0044] (2) The method for correcting the radiation consistency of aircraft multispectral images designed in the present invention introduces vignetting correction into the radiation consistency correction of drone multispectral images, effectively reducing the problem of uneven radiation conditions of a single image caused by the camera lens, and combining the image enhancement algorithm with the SIFT image registration algorithm. In the face of drone multispectral images acquired under different weather conditions and geographical environments, sufficient samples can be provided for the radiation consistency correction model. In addition, the correction model established based on the RANSAC algorithm in the design effectively eliminates gross errors in the sample set and improves the model accuracy. BRIEF DESCRIPTION OF THE DRAWINGS

[0045] Figure 1 It is a flow chart of the method for correcting the radiation consistency of aircraft multi-spectral images designed by the present invention;

[0046] Figure 2 It is a schematic diagram of a multispectral image to be corrected in an embodiment of the design and application of the present invention;

[0047] Figure 3 It is a schematic diagram of a multi-spectral reference image in an embodiment of the design and application of the present invention;

[0048] Figure 4 It is a schematic diagram of a multi-spectral correction image after correction in the design and application embodiment of the present invention. DETAILED DESCRIPTION

[0049] The specific implementation modes of the present invention will be further described in detail below in conjunction with the accompanying drawings.

[0050] The present invention designs a method for correcting the radiation consistency of aircraft multispectral images. Based on the multispectral reference images obtained by the aircraft shooting the target type area at each preset flight altitude, the radiation consistency correction is realized for each multispectral image to be corrected of the target type area shot by the aircraft at each preset flight altitude during the same flight as the multispectral reference image.

[0051] In the application, regarding the multispectral reference image, a calibration plate is placed on the ground in an area of ​​the target type area that is flat, far away from water, and has no shadows. The aircraft performs multispectral photography of the area at each preset flight altitude to obtain the multispectral reference image of the target type area corresponding to each preset flight altitude. For each multispectral reference image, a radial vignetting correction model is applied to correct the vignetting of the multispectral reference image, and each multispectral reference image is updated.

[0052] The correction method includes: for each multispectral image to be corrected, Figure 1 As shown, perform the following steps A to G.

[0053] Step A: Apply the radial vignetting correction model to perform vignetting correction on the multispectral image to be corrected, update the multispectral image to be corrected, and then proceed to step B.

[0054] Step B: Perform histogram contrast stretching processing on the multispectral image to be corrected to obtain a multispectral stretched image to be corrected, and then proceed to step C.

[0055] In actual application, the above step B specifically performs the following steps B1 to B3.

[0056] Step B1. Count the number of pixels corresponding to each different pixel value in the multispectral image to be corrected, and use the ratio of this number to the total number of pixels in the multispectral image to be corrected to form the pixel point proportion corresponding to each different pixel value, and then for each different pixel value, obtain the sum of the pixel point proportions corresponding to each pixel value that is not greater than the pixel value as the cumulative pixel point proportion corresponding to the pixel value, and then obtain the cumulative pixel point proportion corresponding to each different pixel value, and then construct a histogram with the horizontal axis as the pixel value and the vertical axis as the cumulative pixel point proportion, and enter step B2.

[0057] Step B2. Based on the cumulative pixel percentages in the histogram, a cutoff value comparison method is applied to obtain lower pixel values ​​and upper pixel values ​​in different pixel values ​​in the multispectral image to be corrected, and then proceed to step B3.

[0058] In a specific implementation, two different methods are designed for step B2 here to achieve the acquisition of the lower pixel value and the upper pixel value in each different pixel value in the multispectral image to be corrected.

[0059] The first method is to execute step B2-I. Based on the cumulative pixel ratios in the histogram, combined with the preset ratio cutoff value a and b=100%-a, determine whether there is a cumulative pixel ratio equal to a. If so, use the cumulative pixel ratio as the target cumulative pixel ratio, otherwise select the cumulative pixel ratio closest to a as the target cumulative pixel ratio; at the same time, determine whether there is a cumulative pixel ratio equal to b. If so, use the cumulative pixel ratio as the target cumulative pixel ratio, otherwise select the cumulative pixel ratio closest to b as the target cumulative pixel ratio; then obtain the pixel values ​​corresponding to the two target cumulative pixel ratios, use the smaller pixel value of the two pixel values ​​as the lower pixel value and the larger pixel value as the upper pixel value, and then enter step B3.

[0060] In the above step B2-I, when the preset proportion cutoff value a is less than 50%: if the number of cumulative pixel point proportions closest to a is 2, the smaller cumulative pixel point proportion of the two cumulative pixel point proportions is selected as the target cumulative pixel point proportion, and the pixel value corresponding to the target cumulative pixel point proportion is obtained as the lower pixel value; if the number of cumulative pixel point proportions closest to b is 2, the larger cumulative pixel point proportion of the two cumulative pixel point proportions is selected as the target cumulative pixel point proportion, and the pixel value corresponding to the target cumulative pixel point proportion is obtained as the upper pixel value.

[0061] In practical applications, the first method is to set a to 2%, based on the cumulative pixel percentages in the histogram, combined with the preset percentage cutoff values ​​of 2% and 98%, to determine whether there is a cumulative pixel percentage equal to a.

[0062] When the preset proportion cutoff value a is greater than 50%: if the number of cumulative pixel point proportions closest to a is 2, the larger cumulative pixel point proportion of the two cumulative pixel point proportions is selected as the target cumulative pixel point proportion, and the pixel value corresponding to the target cumulative pixel point proportion is obtained as the upper pixel value; if the number of cumulative pixel point proportions closest to b is 2, the smaller cumulative pixel point proportion of the two cumulative pixel point proportions is selected as the target cumulative pixel point proportion, and the pixel value corresponding to the target cumulative pixel point proportion is obtained as the lower pixel value.

[0063] The second method is to execute step B2-II, based on the cumulative pixel percentages sorted from small to large in the histogram, and obtain the number N of the cumulative pixel percentages, combined with the preset percentage threshold c less than 50%, according to Select the dth cumulative pixel point ratio in sequence as the target cumulative pixel point ratio, and select the dth cumulative pixel point ratio from the bottom as the target cumulative pixel point ratio, and obtain the pixel values ​​corresponding to the two target cumulative pixel point ratios respectively, use the smaller pixel value of the two pixel values ​​as the lower pixel value and the larger pixel value as the upper pixel value, and then enter step B3.

[0064] Step B3. Use the lower pixel value and the upper pixel value to replace the maximum pixel value b and the minimum pixel value a in the following formula for contrast stretching of the multispectral image to be corrected, respectively.

[0065]

[0066] Then, the pixel value of each pixel point in the multispectral image to be corrected is updated by the above formula to realize contrast stretching of the multispectral image to be corrected, and the multispectral stretched image to be corrected is obtained, and then step C is entered, wherein f(x, y) represents the pixel value of the pixel point at position (x, y) in the multispectral image to be corrected, g(x, y) represents the pixel value after the pixel value of the pixel point at position (x, y) in the multispectral image to be corrected is updated, and c and d are 0 and 255 respectively.

[0067] Step C: Identify each feature point in the multispectral stretched image to be corrected, and extract the feature vectors corresponding to each feature point, and then enter step D.

[0068] In practical applications, the SIFT algorithm is applied in the above step C, and steps C1 to C5 are specifically performed as follows.

[0069] Step C1. Based on the pixel values ​​of the pixels in the image, a Gaussian function is applied to blur and downsample the multispectral stretched image to be corrected, and an image Gaussian pyramid is constructed, and then step C2 is entered.

[0070] Step C2. In the obtained Gaussian difference pyramid space, for each pixel point, the pixel point is compared with its 8 adjacent pixel points and 2×9=18 points in the adjacent upper and lower layers, a total of 26 points. If the pixel point is the maximum or minimum, it is used as a local key point; then each local key point is obtained, and then step C3 is entered.

[0071] Step C3. Based on each local key point, the position and scale of each local key point are accurately determined by fitting a three-dimensional quadratic function. At the same time, the edge response points are detected according to the Harris algorithm and removed from all local key points to obtain the feature points in the multispectral stretched image to be corrected, and then enter step C4.

[0072] Step C4. For each feature point in the multispectral stretched image to be corrected, obtain the gradient between the feature point and its adjacent pixel points in each direction, and select the direction of the pixel point corresponding to the maximum gradient as the direction corresponding to the feature point; then obtain the directions corresponding to each feature point, and then enter step C5.

[0073] Step C5. For each feature point in the multispectral stretched image to be corrected, a feature vector corresponding to the feature point is constructed using the direction corresponding to the feature point, the position of the feature point in the image, and the gradient between the feature point and its adjacent pixel points in each direction; and then the feature vectors corresponding to each feature point in the multispectral stretched image to be corrected are obtained.

[0074] Step D. Based on the feature points in the multispectral reference image with the same flight altitude as the multispectral image to be corrected obtained in advance according to the method of step C, and the feature vectors corresponding to the feature points, by comparing the Euclidean distances between the feature points, the feature point groups that match each other in the multispectral stretched image to be corrected and the multispectral reference image are obtained, and then the corresponding feature point position groups are obtained, and then step E is entered.

[0075] In the actual implementation design, in the above step D, based on the feature points in the multispectral reference image with the same flight altitude as the multispectral image to be corrected, which are obtained in advance according to the method of step C, and the feature vectors corresponding to each feature point, for each feature point in the multispectral stretched image to be corrected, the feature point is used as a candidate feature point, and the following steps D1 to D2 are performed to realize the search of feature points matching the candidate feature points in the multispectral reference image, obtain the feature point groups that match each other in the multispectral stretched image to be corrected and the multispectral reference image, and the corresponding feature point position groups, and then enter step E.

[0076] Step D1. Obtain the Euclidean distance between each feature point based on the feature vector corresponding to the feature point, obtain the Euclidean distance between the feature points to be selected and each feature point in the multispectral reference image, and obtain the ratio between the closest distance and the second closest distance, and then enter step D2.

[0077] Step D2. Determine whether the ratio is less than a preset distance ratio threshold. If so, determine that the feature point to be selected matches the feature point in the multispectral reference image corresponding to the nearest distance to each other to form a group of feature points, and obtain the corresponding position of the group of feature points; otherwise, determine that there is no feature point in the multispectral reference image that matches the feature point to be selected.

[0078] Step E. For each feature point position group, respectively, obtain the pixel value of the corresponding position in the corresponding feature point position group in the multispectral reference image with the same flight altitude as the multispectral image to be corrected as the reference pixel value, and obtain the pixel value of the corresponding position in the corresponding feature point position group in the multispectral image to be corrected as the pixel value to be corrected, and the pixel value to be corrected and the reference pixel value are combined to form a pixel fitting sample group corresponding to the feature point position group; then obtain the pixel fitting sample groups corresponding to each feature point position group, and then enter step F.

[0079] Step F. Based on the correspondence between the pixel values ​​to be corrected and the reference pixel values ​​in each pixel fitting sample group, data fitting is performed for the preset model according to the random sampling consistency algorithm RANSAC to construct an image radiation consistency correction model for the flight altitude of the aircraft corresponding to the multispectral image to be corrected, and then enter step G.

[0080] Step G: Apply the image radiation consistency correction model to correct and update the pixel values ​​at each position in the multispectral image to be corrected, obtain the multispectral correction image corresponding to the multispectral image to be corrected, and then enter step H.

[0081] Step H: performing bilateral filtering on the multi-spectral correction image to update the multi-spectral correction image.

[0082] The aircraft multispectral image radiation consistency correction method designed by the present invention is applied in practice. In a citrus plantation with water, a DJI Phantom 4-Multispectral UAV (P4M) is used to fly in a sunny and cloudless weather around 14:00 noon at an altitude of 120m. The lateral and directional overlap rates of the captured images are set to 60% and 80% respectively, and the camera AE is set to unlocked state. According to the designed consistency correction method, the following steps are specifically performed:

[0083] Step A: Apply the radial vignetting correction model to Figure 2 The multispectral image acquired by P4M is corrected for vignetting and the multispectral image is updated. The radial vignetting correction model is defined as follows:

[0084] L(x,y)=V(x,y)×P(x,y)

[0085] Where (x, y) is the coordinate of the pixel in the multispectral image, L(x, y) is the pixel value of the (x, y) coordinate after vignetting correction of the multispectral image, and V(x, y) is the gain function of vignetting compensation defined as follows:

[0086] V(x,y)=1+k 0 r+k 1 r 2+....+k 4 r 5 +k 5 r 6

[0087] Among them, k 0 ,k 1 ,....,k 4 ,k 5 are optical parameters. These hardware-related parameters are given by the manufacturer and can be read in the image metadata information. r is the pixel distance from the (x, y) coordinate pixel point in the multispectral image to the compensation center. The calculation formula is as follows:

[0088]

[0089] CenterX and CenterY are the image centers for vignetting compensation.

[0090] Step B: Perform histogram contrast stretching processing on the multispectral image to be corrected to obtain the multispectral stretched image to be corrected. The histogram contrast stretching algorithm is defined as follows:

[0091]

[0092] Where: f(x,y) represents the pixel value of the pixel at position (x,y) in the multispectral image to be corrected, g(x,y) represents the pixel value after updating the pixel value of the pixel at position (x,y) in the multispectral image to be corrected, a and b are the minimum and maximum values ​​of the original image pixel values, respectively. The preset percentage cutoff value a is set to 2%, so the c and d values ​​are the pixel values ​​corresponding to 2% and 98% in the cumulative histogram.

[0093] Step C: using the SIFT algorithm to identify each feature point in the multi-spectral stretched image to be corrected, and extracting the feature vectors corresponding to each feature point.

[0094] Step D: Based on Figure 3 The multispectral reference image shown is combined with each feature point in the multispectral reference image obtained in advance according to the method in step C at the same flight altitude as the multispectral image to be corrected, to obtain each feature point group in the multispectral stretched image to be corrected that matches each other in the multispectral reference image; wherein, the threshold value of the ratio of the Euclidean distance of the nearest neighbor feature point to the Euclidean distance of the next nearest neighbor feature point between a feature point of the multispectral stretched image to be corrected and the multispectral reference image is set to 0.5, that is, less than 0.5, it is considered that the two feature points match each other, constitute a set of feature point pairs, and obtain the corresponding position of the set of feature points; otherwise, it is determined that there is no feature point in the multispectral reference image that matches the feature point in the multispectral stretched image to be corrected.

[0095] Step E. For each feature point position group obtained in step D, extract the pixel values ​​of the corresponding positions in the feature point position group of the reference image as reference pixel values, and extract the pixel values ​​of the corresponding positions in the corresponding feature point position group in the image to be corrected as pixel values ​​to be corrected, and form a pixel fitting sample group corresponding to the feature point position group by combining the pixel values ​​to be corrected and the reference pixel values; and then obtain the pixel fitting sample groups corresponding to each feature point position group.

[0096] Step F: Based on the correspondence between the pixel values ​​to be corrected and the reference pixel values ​​in each pixel fitting sample group, linear regression is performed according to the RANSAC algorithm to establish a radiation consistency correction model between the image to be corrected and the reference image.

[0097] In this embodiment, the step F specifically includes the following steps:

[0098] F1 randomly selects a certain number of samples as internal points;

[0099] F2 calculates the model suitable for the inliers;

[0100] F3 uses the model to test all other points and puts the points that fall within the given range into the inner point set;

[0101] F4 records the number of internal points;

[0102] F5 Repeat the above steps multiple times;

[0103] F6 takes the model with the largest number of inliers as the final regression model, where the model function is as follows:

[0104] y i =kx i +bi=1,2,3,...,n

[0105] Among them, x i is the pixel value to be corrected in the i-th pixel fitting sample group, y i is the reference pixel value in the i-th pixel fitting sample group, and n is the total number of pixel fitting sample groups.

[0106] Step G: Apply the image radiation consistency correction model to correct and update the pixel value at each position in the multispectral image to be corrected, and obtain the multispectral correction image corresponding to the multispectral image to be corrected.

[0107] Step H: Perform bilateral filtering on the multispectral correction image to update the multispectral correction image, such as Figure 4 As shown. The bilateral filtering model is defined as follows:

[0108]

[0109] Among them, (i,j) is the image pixel coordinate, g(i,j) is the pixel value of the output pixel at position (i,j), f(k,l) is the pixel value of the pixel in the input filter at position (i,j) of the original image, (k,l) is the filter size, and w(i,j,k,l) ​​is the weighting coefficient determined by the product of the spatial kernel and the range kernel:

[0110] w(i,j,k,l)=d(i,j,k,l)×r(i,j,k,l)

[0111] The spatial kernel expression is as follows:

[0112]

[0113] The range kernel expression is as follows:

[0114]

[0115] The filter size used in this embodiment is 3×3, and the standard deviation of the bilateral filter space domain σ d =10, grayscale standard deviation σ r =30 (grayscale range [0~255]).

[0116] The radiation consistency correction method for aircraft multispectral images designed by the above technical scheme first performs vignetting correction; then performs histogram contrast stretching on the image; then uses the SIFT (Scale invariant feature transform) operator to perform feature matching on the reference image and the image after histogram contrast stretching to obtain matching point pairs; then based on the pixel values ​​of the matching point pairs, linear regression is performed according to the RANSAC (Random sample consensus) algorithm to establish a linear correction model; finally, the linear correction model is applied to perform radiation consistency correction on the image to be corrected, and bilateral filtering is performed on the corrected image to remove noise points; the design scheme overcomes the disadvantage that the traditional correction method cannot effectively obtain sufficient same-name points in low-illuminance images and cannot carry out radiation information correction between images; at the same time, the modeling method has strong robustness and effectively resists the interference of noise points in the image, which can make the radiation conditions of thousands of drone images in a sortie tend to be consistent, eliminating the image distortion problem caused by strong reflection of ground objects and drastic changes in light intensity; thereby providing solid data support for subsequent remote sensing quantitative research.

[0117] In practical applications, vignetting correction is introduced into the radiation consistency correction of UAV multispectral images, which effectively reduces the problem of uneven radiation conditions in a single image caused by the camera lens, and combines the image enhancement algorithm with the SIFT image registration algorithm. In the face of UAV multispectral images acquired under different weather conditions and geographical environments, they can provide sufficient samples for the radiation consistency correction model. In addition, the correction model established based on the RANSAC algorithm in the design effectively eliminates gross errors in the sample set and improves the model accuracy.

[0118] The embodiments of the present invention are described in detail above with reference to the accompanying drawings, but the present invention is not limited to the above embodiments, and various changes can be made within the knowledge scope of ordinary technicians in this field without departing from the purpose of the present invention.

Claims

1. A method for correcting the radiometric consistency of aircraft multispectral images. Features: Based on the multispectral reference images obtained by the aircraft photographing the target type area at each preset flight altitude, radiation consistency correction is performed for each multispectral image to be corrected of the target type area photographed by the aircraft at each preset flight altitude during the same flight as the multispectral reference image. The correction method includes performing the following steps A to G for each multispectral image to be corrected; Step A. Apply the radial vignetting correction model to perform vignetting correction on the multispectral image to be corrected, update the multispectral image to be corrected, and then proceed to step B; The radial vignetting correction model in step A is defined as follows: L(x,y)=V(x,y)×P(x,y) Where (x, y) is the coordinate of the pixel in the multispectral image, L(x, y) is the pixel value of the (x, y) coordinate after vignetting correction of the multispectral image, and V(x, y) is the gain function of vignetting compensation defined as follows: V(x,y)=1+k 0 r+k 1 r 2 +....+k 4 r 5 +k 5 r 6 Among them, k 0 ,k 1 ,....,k 4 ,k 5 are optical parameters. These hardware-related parameters are given by the manufacturer and can be read in the image metadata information. r is the pixel distance from the (x, y) coordinate pixel point in the multispectral image to the compensation center. The calculation formula is as follows: CenterX and CenterY are the image centers for vignetting compensation; Step B. Performing histogram contrast stretching processing on the multispectral image to be corrected to obtain the multispectral stretched image to be corrected, and then proceeding to step C; Step C. using the SIFT algorithm to identify each feature point in the multispectral stretched image to be corrected, and extracting the feature vectors corresponding to each feature point, and then proceeding to step D; Step D. Based on the feature points in the multispectral reference image with the same flight altitude as the multispectral image to be corrected obtained in advance according to the method in step C, and the feature vectors corresponding to the feature points, by comparing the Euclidean distances between the feature points, obtain the feature point groups that match each other in the multispectral stretched image to be corrected and the multispectral reference image, and then obtain the corresponding feature point position groups, and then enter step E; Step E. For each feature point position group, respectively, obtain the pixel value of the corresponding position in the feature point position group in the multispectral reference image with the same flight altitude as the multispectral image to be corrected, as the reference pixel value, and obtain the pixel value of the corresponding position in the feature point position group in the multispectral image to be corrected, as the pixel value to be corrected, and form a pixel fitting sample group corresponding to the feature point position group by combining the pixel value to be corrected and the reference pixel value; then obtain the pixel fitting sample groups corresponding to each feature point position group, and then enter step F; Step F. Based on the correspondence between the pixel values ​​to be corrected and the reference pixel values ​​in each pixel fitting sample group, data fitting is performed for a preset model according to the random sampling consistency algorithm RANSAC to form an image radiation consistency correction model corresponding to the flight altitude of the aircraft and the multispectral image to be corrected, and then proceeding to step G; Step G. Apply the image radiation consistency correction model to correct and update the pixel values ​​at each position in the multispectral image to be corrected, obtain the multispectral correction image corresponding to the multispectral image to be corrected, and then enter step H; Step H: performing bilateral filtering on the multispectral correction image to update the multispectral correction image.

2. According to claim 1, a method for correcting the radiation consistency of aircraft multispectral images, Features: A calibration plate is placed on the ground in an area of ​​the target type area that is flat, far away from water and has no shadows. The aircraft takes multispectral photos of the area at each preset flight altitude to obtain multispectral reference images of the target type area corresponding to each preset flight altitude. For each multispectral reference image, a radial vignetting correction model is applied to correct the vignetting of the multispectral reference image, and each multispectral reference image is updated.

3. According to claim 1, a method for correcting the radiation consistency of aircraft multispectral images, Features: The step B includes steps B1 to B3 as follows: Step B1. Count the number of pixels corresponding to each different pixel value in the multispectral image to be corrected, and form the pixel ratio corresponding to each different pixel value through the ratio of the number to the total number of pixels in the multispectral image to be corrected, and then obtain the sum of the pixel ratios corresponding to each pixel value not greater than the pixel value for each different pixel value as the cumulative pixel ratio corresponding to the pixel value, and then obtain the cumulative pixel ratio corresponding to each different pixel value, and then construct a histogram with the horizontal axis as the pixel value and the vertical axis as the cumulative pixel ratio, and enter step B2; Step B2. Based on the cumulative pixel percentages in the histogram, a cutoff value comparison method is applied to obtain the lower pixel value and the upper pixel value of each different pixel value in the multispectral image to be corrected, and then proceed to step B3; Step B3. Use the lower pixel value and the upper pixel value to replace the maximum pixel value and the minimum pixel value in the formula involved in contrast stretching of the multispectral image to be corrected, and then perform contrast stretching on the multispectral image to be corrected to obtain the multispectral stretched image to be corrected, and then enter step C.

4. According to claim 3, a method for correcting the radiation consistency of aircraft multispectral images, Features: The step B2 is the following step B2-I; Step B2-I. Based on the cumulative pixel percentages in the histogram, combined with the preset percentage cutoff value a and b=100%-a, determine whether there is a cumulative pixel percentage equal to a, if yes, use the cumulative pixel percentage as the target cumulative pixel percentage, otherwise select the cumulative pixel percentage closest to a as the target cumulative pixel percentage; at the same time, determine whether there is a cumulative pixel percentage equal to b, if yes, use the cumulative pixel percentage as the target cumulative pixel percentage, otherwise select the cumulative pixel percentage closest to b as the target cumulative pixel percentage; Then, the pixel values ​​corresponding to the two target cumulative pixel point proportions are obtained, and the smaller pixel value of the two pixel values ​​is used as the lower pixel value and the larger pixel value is used as the upper pixel value, and then step B3 is entered.

5. According to claim 4, a method for correcting the radiation consistency of aircraft multispectral images, Features: In step B2-I, when the preset percentage cutoff value a is less than 50%: If the number of cumulative pixel ratios closest to a is 2, the smaller cumulative pixel ratio of the two cumulative pixel ratios is selected as the target cumulative pixel ratio, and the pixel value corresponding to the target cumulative pixel ratio is obtained as the next pixel value; If the number of cumulative pixel ratios closest to b is 2, the larger cumulative pixel ratio of the two cumulative pixel ratios is selected as the target cumulative pixel ratio, and the pixel value corresponding to the target cumulative pixel ratio is obtained as the upper pixel value; When the preset percentage cutoff value a is greater than 50%: If the number of cumulative pixel ratios closest to a is 2, the larger cumulative pixel ratio of the two cumulative pixel ratios is selected as the target cumulative pixel ratio, and the pixel value corresponding to the target cumulative pixel ratio is obtained as the upper pixel value; If the number of cumulative pixel point ratios closest to b is 2, the smaller cumulative pixel point ratio of the two cumulative pixel point ratios is selected as the target cumulative pixel point ratio, and the pixel value corresponding to the target cumulative pixel point ratio is obtained as the next pixel value.

6. According to claim 3, a method for correcting the radiation consistency of aircraft multispectral images, Features: The step B2 is the following step B2-II; Step B2-II. Based on the cumulative pixel percentages in the histogram sorted from small to large, the number N of cumulative pixel percentages is obtained, combined with the preset percentage threshold c of less than 50%, according to Select the dth cumulative pixel point ratio in sequence as the target cumulative pixel point ratio, and select the dth cumulative pixel point ratio from the bottom as the target cumulative pixel point ratio, and obtain the pixel values ​​corresponding to the two target cumulative pixel point ratios respectively, use the smaller pixel value of the two pixel values ​​as the lower pixel value and the larger pixel value as the upper pixel value, and then enter step B3.

7. According to claim 3, a method for correcting the radiation consistency of aircraft multispectral images, Features: In step B3, the lower pixel value and the upper pixel value are used to replace the maximum pixel value b and the minimum pixel value a in the following formula for contrast stretching of the multispectral image to be corrected, respectively. Then, the pixel value of each pixel point in the multispectral image to be corrected is updated by the above formula to realize contrast stretching of the multispectral image to be corrected, and obtain the multispectral stretched image to be corrected, wherein f(x, y) represents the pixel value of the pixel point at position (x, y) in the multispectral image to be corrected, g(x, y) represents the pixel value after the pixel value of the pixel point at position (x, y) in the multispectral image to be corrected is updated, and c and d are 0 and 255 respectively.

8. According to claim 1, a method for correcting the radiation consistency of aircraft multispectral images, Features: In step C, the SIFT algorithm is applied, and steps C1 to C5 are performed as follows: Step C1. Based on the pixel values ​​of the pixels in the image, a Gaussian function is applied to blur and downsample the multispectral stretched image to be corrected, and an image Gaussian pyramid is constructed, and then step C2 is entered; Step C2. In the obtained Gaussian difference pyramid space, for each pixel point, the pixel point is compared with its 8 adjacent pixel points and 2×9=18 points in the adjacent upper and lower layers, a total of 26 points. If the pixel point is the maximum or minimum, it is used as a local key point; then each local key point is obtained, and then step C3 is entered; Step C3. Based on each local key point, the position and scale of each local key point are accurately determined by fitting a three-dimensional quadratic function, and the edge response point is detected according to the Harris algorithm, and removed from all local key points, thereby obtaining the feature points in the multispectral stretched image to be corrected, and then entering step C4; Step C4. For each feature point in the multispectral stretched image to be corrected, obtain the gradient between the feature point and its adjacent pixel points in each direction, and select the direction of the pixel point corresponding to the maximum gradient as the direction corresponding to the feature point; then obtain the directions corresponding to each feature point, and then enter step C5; Step C5. For each feature point in the multispectral stretched image to be corrected, a feature vector corresponding to the feature point is constructed using the direction corresponding to the feature point, the position of the feature point in the image, and the gradient between the feature point and its adjacent pixel points in each direction; and then the feature vectors corresponding to each feature point in the multispectral stretched image to be corrected are obtained.

9. According to claim 1, a method for correcting the radiation consistency of aircraft multispectral images, Features: In the step D, based on each feature point in the multispectral reference image with the same flight altitude as the multispectral image to be corrected obtained in advance according to the method of step C, and the feature vectors corresponding to each feature point, each feature point in the multispectral stretched image to be corrected is used as a candidate feature point, and the following steps D1 to D2 are performed to search for feature points matching the candidate feature points in the multispectral reference image, obtain each feature point group matching each other in the multispectral stretched image to be corrected and the multispectral reference image, and each corresponding feature point position group, and then enter step E; Step D1. Obtain the Euclidean distance between each feature point in the multispectral reference image by obtaining the Euclidean distance between each feature point based on the feature vector corresponding to the feature point, and obtain the ratio between the closest distance and the second closest distance, and then proceed to step D2; Step D2. Determine whether the ratio is less than a preset distance ratio threshold. If so, determine that the feature point to be selected matches the feature point in the multispectral reference image corresponding to the nearest distance to each other to form a group of feature points, and obtain the corresponding position of the group of feature points; otherwise, determine that there is no feature point in the multispectral reference image that matches the feature point to be selected.

Images