CN108961162A - A kind of unmanned plane forest zone Aerial Images joining method and system - Google Patents

Abstract

The invention discloses a kind of forest zone Aerial Images joining method and system based on unmanned plane, comprising the following steps: pass through the visible light collection equipment of UAV flight and shoot forest zone image；To the forest zone Aerial Images that the visible light collection equipment obtains, image preprocessing is carried out using color histogram equalization algorithm, picture contrast is improved, enhances image detail information, increase feature point extraction number；Extreme point detection and extraction are then carried out by the contrast threshold screening method of optimization, the extraction of unnecessary characteristic point is reduced, improves detection efficiency, while replacing the ratio between Euclidean distance to carry out feature using the angle ratio of dot product inverse cosine function and slightly matching, computation complexity is reduced, matching efficiency is improved；Then it on the basis of using RANSAC algorithm removal Mismatching point, introduces L-M non-linear iterative and refines inter-image transformations matrix, carry out fine match, improve image registration accuracy；Image mosaic is finally realized using weighted average blending algorithm.

Description

A method and system for splicing aerial images of unmanned aerial vehicle forest areas

technical field

The invention relates to the field of stitching and processing of aerial images, in particular to a method and system for stitching images of aerial images of unmanned forest areas.

Background technique

As an emerging low-altitude remote sensing technology in recent years, UAV aerial photography has the advantages of strong timeliness, low cost, and high safety. It has been widely used in forestry surveys such as forestry pest monitoring, forest fire monitoring, and forest stand analysis.

However, due to the limitations of flight height, digital camera focal length, high spatial resolution and other factors when the UAV remote sensing platform acquires images, the coverage of a single UAV image is limited. Therefore, it is necessary to stitch the aerial images of UAVs to provide high-resolution, wide-field aerial images of forest areas for the analysis and processing of remote sensing images for forestry investigation and research.

When the traditional aerial image stitching algorithm processes aerial images in forest areas, there are problems such as too few feature points to be extracted, complex matching methods, and low registration accuracy. In order to solve this problem, this patent proposes an IE-SIFT ( ImageEnhanced Scale-Invariant Feature Transform) image stitching algorithm. The algorithm first introduces global histogram equalization in the image preprocessing stage to improve image contrast, which enhances image detail information and increases the number of feature points to be extracted; then uses the optimized contrast threshold screening method to detect and extract extreme points to reduce unnecessary features Point extraction improves detection efficiency. At the same time, the angle ratio of the vector dot product arccosine function is used instead of the ratio of Euclidean distance for feature rough matching, which reduces computational complexity and improves matching efficiency; Firstly, the L-M non-linear iterative algorithm is introduced to refine the transformation matrix between images, for fine matching, and to improve the accuracy of image registration; finally, the weighted average fusion algorithm is used to realize image mosaic.

Contents of the invention

In view of this, the object of the present invention is to propose a method and system for splicing aerial images of UAV forest areas to realize fast and effective splicing of aerial images of UAV forest areas.

Based on the above-mentioned purpose, the method and system for splicing aerial images of unmanned aerial vehicles in forest areas provided by the present invention include steps:

Take images of the forest area over the target area through the drone equipped with visible light acquisition equipment;

Using a global histogram equalization algorithm to preprocess the aerial image to obtain an enhanced grayscale aerial image;

Through the optimized contrast threshold screening method, the enhanced grayscale aerial image is subjected to feature point detection and extraction to reduce the extraction of unnecessary feature points, and an improved feature point matching method is used to compare the reference image and the image to be registered The feature points are matched, and the schematic diagram of feature point matching is obtained;

Use the RANSAC algorithm to remove mismatching points and calculate the initial value of the transformation matrix, use the L-M nonlinear iterative algorithm to refine the calculated initial transformation matrix, and continuously correct the values of the parameters of the transformation matrix through iterative calculations to obtain the final accurate transformation Matrix, using the final transformation matrix to obtain the image to be fused after performing image transformation on the original image of the image to be registered;

performing image fusion on the original image of the reference image and the image to be fused by using a weighted average fusion algorithm to obtain a final fused image.

As an embodiment, the overlapping degree of the captured images is 40%-60%.

As an embodiment, the processing process of the global histogram equalization algorithm is as follows: normalize the gray level r of the aerial image _{gzh(x, y)} to the interval [0,1], when r=0, it is black, When r=1, it is white; _{gzh(x, y)} gray scale range is [0, L-1], the total number of pixels is n, then the number of pixels with gray scale r _k is _nk ; The transformation corresponding to global histogram equalization is shown in formula (1);

P _r (r _j ) is the probability distribution function of g _{zh(x, y)} , and S _k is when the gray level k of the aerial image g _{zh(x, y)} is taken from 0 to (L-1), Summing the image _{gzh(x, y)} probability distribution function P _r (r _j ), L is the total number of gray levels of the aerial image, T (r) represents a transformation function; use the transformation formula to change the aerial image The grayscale of the enhanced grayscale aerial image is obtained.

As an embodiment, the optimized contrast threshold screening method is as follows: first, the points in the image whose pixel value is lower than the contrast threshold are eliminated, and the method for selecting the threshold is: the correct registration of the threshold size and the feature points is obtained through experimental verification. ratio, and select the optimal contrast threshold.

As an embodiment, the processing procedure of the improved feature point matching algorithm is as follows: perform a dot product operation on a certain feature point in the reference image and all the feature points in the image to be stitched, and then perform an inverse cosine transform on the dot product result to store to the specified array, and then find the minimum angle θ ₁ and the second minimum angle θ ₂ from the array, if the ratio of the two is less than the specified threshold M=0.6, then it is considered that the feature point e _l corresponding to the minimum angle is the same as that in the reference image match the feature points e _r ; (e _l , e _r ) is left as a matching pair of feature points, and the matching formula is shown in formula (2):

In this way, initial matching can be performed on the extracted feature points to obtain initial matching point pairs.

As an embodiment, the use of RANSAC algorithm to remove mismatching points is using formula (3):

Remove the mismatching points to obtain corresponding matching points; in formula (3), H is a transformation matrix, and the transformation matrix H is a 3×3 matrix containing 8 unknown parameters (m ₁₁ ~m ₃₂ ); (x _i , y _i ), (x′ _i , y′ _i ) are corresponding matching points A _l and A _r of the reference image and the image to be stitched, respectively.

As an embodiment, the L-M nonlinear iterative algorithm is used to refine the calculated initial transformation matrix, and the value of each parameter of the transformation matrix is continuously corrected by iterative calculation to obtain the final accurate transformation matrix, including steps:

S601, set n=0; μ=0.01; set a threshold τ; the maximum number of iterations is N; the initial transformation matrix is m ₀ ;

S602, calculating the value of F(m _n );

S603, judge whether F(m _n ) is greater than τ, if so, go to S604, otherwise end the process;

S604, judge n! Whether it is greater than 0, if so, go to S605, otherwise go to S606;

S605, judge whether F(m _n ) is smaller than F(m _n-1 ), if so, go to S606, otherwise, set m _n =m _n-1 , μ=10×μ and go to S607;

S606, set μ=μ/10 and calculate the value of J(m _n ), go to S607;

S607, judge whether n is greater than N, if so, end the process, otherwise, according to the formula:

m _n+1 ＝m _n －[J(m _n ) ^τ +J(m _n )μI] ^－1 J(m _n ) ^τ f(m _n ), where J(m) is the Jacobian matrix, let n =n+1, go to S602, namely calculate the value of F(m _n );

Use the value of F(m _n ) into (4)-Equation (6)

In the formula (m ₁₁ ~m ₃₂ ) are the components of the transformation matrix H, and in the formula (6), F(m) is the sum of the distances between all corresponding matching feature points; the optimal solution of each component is calculated as The final transformation moment H _b is obtained.

A method and system for mosaicing aerial images of unmanned aerial vehicle forest areas, comprising:

The image input module is used to obtain the forest area aerial image over the target area;

An image enhancement module, configured to use a global histogram equalization algorithm to preprocess the aerial image to obtain an enhanced grayscale aerial image;

The feature point extraction module is used to use the optimized contrast threshold screening method to perform scale space extremum detection on the enhanced gray-scale aerial image, accurately locate feature points, determine the main direction of feature points, and generate feature point descriptors; use the The feature point descriptor is used to determine the feature point, and the reference image and the image to be registered of the feature point are obtained;

The feature point matching module is used to perform feature point matching on the feature points of the reference image and the image to be registered using the ratio of the arccosine function of the unit vector dot product to obtain a schematic diagram of feature point matching;

Remove the mismatching point module, which is used to make the RANSAC algorithm remove the mismatching points in the feature point matching schematic diagram, obtain the matching result schematic diagram after removing the mismatching points and calculate and obtain the initial transformation matrix;

The refining transformation matrix module is used to refine the initial transformation matrix using the L-M nonlinear iterative algorithm, so that the sum of the distances between all corresponding matching feature points in the matching result schematic diagram after removing the mismatching points is the smallest, and obtains the final transformation matrix; Using the final transformation matrix to obtain the image to be fused after performing image transformation on the original image of the image to be registered;

An image fusion module, configured to use a weighted average fusion algorithm to perform image fusion on the original image of the reference image and the image to be fused to obtain a final fused image.

Description of drawings

In order to more clearly illustrate the technical solutions in the embodiments of the present invention or the prior art, the following will briefly introduce the drawings that need to be used in the description of the embodiments or the prior art. Obviously, the accompanying drawings in the following description are only These are some embodiments of the present invention. Those skilled in the art can also obtain other drawings based on these drawings without creative work.

Fig. 1 is a kind of flow chart of the splicing method of unmanned aerial vehicle forest area aerial photography image of the embodiment of the present invention;

Fig. 2a is the original aerial image of the embodiment of the present invention;

Figure 2b is a grayscale aerial image of an embodiment of the present invention;

FIG. 2c is an enhanced grayscale aerial image of an embodiment of the present invention;

Fig. 3a is a schematic diagram of a feature point extraction result according to an embodiment of the present invention;

Fig. 3b is a schematic diagram of another feature point extraction result according to the embodiment of the present invention;

FIG. 4 is a schematic diagram of a feature point matching result according to an embodiment of the present invention;

FIG. 5 is a schematic diagram of a matching result after removing mismatch points according to an embodiment of the present invention;

Fig. 6 is a flow chart of L-M algorithm refining transformation matrix according to the embodiment of the present invention;

Fig. 7a is the original image of the image to be registered according to the embodiment of the present invention;

Fig. 7b is an image to be fused according to an embodiment of the present invention;

Fig. 8 is the final fused image of the embodiment of the present invention;

Fig. 9 is a structural diagram of an aerial image mosaic system of an unmanned aerial vehicle forest area according to an embodiment of the present invention.

Detailed ways

In order to make the object, technical solution and advantages of the present invention clearer, the present invention will be described in further detail below in conjunction with specific embodiments and with reference to the accompanying drawings.

As shown in Figure 1, it is a flow chart of a method for splicing aerial photographs of unmanned aerial vehicles in forest areas according to an embodiment of the present invention, and the specific implementation steps include:

In step 101, images are obtained, and the drone will take aerial images of forest areas during its flight over the target area.

Step 102, histogram equalization image enhancement, because the aerial image in some areas has a single color and the change is not obvious, resulting in the conversion of the aerial image into a grayscale image, the contrast is too low, the image details are blurred, and it is difficult to effectively use the SIFT algorithm to extract features point. Therefore, after the image processing terminal receives the aerial image, the global histogram equalization algorithm is used to preprocess the aerial image, which can improve image contrast, increase image detail information, and thereby increase the number of extracted feature points.

The specific processing process is as follows: the gray level r of the aerial image _gzh(x,y) is normalized to the interval [0,1], when r=0 it is black, and when r=1 it is white. g _zh(x,y) gray scale range is

[0,L-1], the total number of pixels is n, then the number of pixels with gray level r _k is n _k . The transformation corresponding to the global histogram equalization is shown in formula (1).

In formula (1), P _r (r _j ) is the probability distribution function of g _{zh(x, y)} , and S _k is the gray level k of the aerial image g _{zh(x, y)} from 0 to (L -1), the probability distribution function P _r (r _j ) of the image _gzh(x,y) is summed, L is the total number of gray levels of the aerial image, and T(r) represents the transformation function. Using the transformation formula to change the grayscale of the aerial image to obtain an enhanced grayscale aerial image. The contrast effect after image enhancement is shown in Figure 2a-2c. The contrast of the enhanced grayscale aerial image is significantly improved, and the image detail information is more abundant. 2a is the original aerial image of the embodiment of the present invention, FIG. 2b is the grayscale aerial image of the embodiment of the present invention, and FIG. 2c is the enhanced grayscale aerial image of the embodiment of the present invention.

Step 103, using the SIFT algorithm to extract feature points. When the UAV is shooting the aerial image of the target area, due to the interference of illumination changes and other environmental factors, the input image sequence may be rotated. The perspective and scale of the aerial image Changes may also occur. In order to overcome the above interference, the embodiment of the present invention uses the SIFT algorithm to perform image registration. The image registration is the core step of image stitching, and the accuracy of registration will directly affect the final effect of image stitching. The image registration is to calculate the spatial transformation parameters between the reference image and the image to be registered, determine the spatial transformation model, and realize the transformation of the coordinate system between the two images. The SIFT algorithm has good scale invariance and robustness to viewing angle changes by extracting feature points in the scale space.

The extracted feature points include:

Scale space extreme value detection, the scale space L(x,y,σ) of a two-dimensional image is defined as the convolution of the scale-variable Gaussian function G(x,y,σ) and the original image I(x,y), As shown in formulas (2) and (3).

L(x,y,σ)=G(x,y,σ)*I(x,y) (3)

In formula (3), * is the convolution operation, (x, y) is the position of the image pixel, σ is the scale space factor, the larger the value of σ, the blurrier the image, and the smaller the value of σ, the clearer the image. In order to effectively detect relatively stable feature points in the scale space, the Gaussian difference scale space D(x, y, σ) is generated by convolution of Gaussian difference kernels of different scales with the image. As shown in formula (4).

D(x,y,σ)=(G(x,y,kσ)-G(x,y,σ))*I(x,y)=L(x,y,kσ)-L(x,y ,σ)

(4)

In formula (4), k is a constant multiplicative factor, and its value is related to the number S of pictures per layer in the scale space, that is, k=2 ^1/S . The enhanced grayscale aerial image is detected, and if the grayscale value of a certain pixel reaches an extreme value in 26 neighborhoods around and up and down, record the position and scale of this point as a candidate feature point. This patent combines the contrast threshold screening method with the features and characteristics of the forest area image to detect the extreme points, reduce the extraction of invalid feature points in the image, and thus improve the efficiency of the algorithm for detecting and extracting the feature points of this type of image. The principle of threshold selection is to ensure the correct registration in the forest image while minimizing feature point extraction. The relationship between the threshold size and the correct registration rate of feature points is obtained through experimental verification, so as to select the optimal contrast threshold. The main basis of this method is: there are a large number of similar structures in the forest image, and the pixel value changes between the forest scenes are small, but there are obvious contrast and pixel value changes between forest trees and non-forest trees in the image. Therefore, the method of increasing the contrast threshold can be used to extract the main scene in the image content, thereby extracting feature points that are more likely to be useful for image registration, and reducing the detection of feature points that do not contribute to registration.

Precisely locate feature points, accurately determine the position and scale of feature points by fitting three-dimensional quadratic functions, and filter low-contrast feature points and unstable edge response points to enhance matching stability and improve noise resistance.

The main direction of the feature point is determined. In order to make the descriptor invariant to rotation, a direction is assigned to each feature point by using the local features of the grayscale aerial image. Using the characteristics of the gradient and direction distribution of the neighborhood pixels of the feature point, the modulus m(x, y) and the direction θ(x, y) of the gradient are obtained, as shown in formulas (5) and (6).

To generate a feature point descriptor, in order to ensure the invariance of feature vector rotation, first rotate the coordinate axis to the main direction of the feature point, use a Gaussian circular window to Gaussian weight the gradient modulus, and use the feature point as the center of the 16×16 window, Divide this window into 4×4 sub-regions, each region counts 8 directions through the gradient histogram, and finally forms a 128-dimensional SIFT feature vector that is not affected by scale, rotation and illumination, which is the feature point descriptor.

Feature points are determined using the feature point descriptor.

As shown in Figure 3a and Figure 3b, it is the feature point extraction result. FIG. 3 a is a schematic diagram of a feature point extraction result according to an embodiment of the present invention, and FIG. 3 b is a schematic diagram of another feature point extraction result according to an embodiment of the present invention.

Step 104, feature point matching. In order to simplify the matching process and improve the feature point matching speed, the embodiment of the present invention uses the ratio of the arccosine function of the unit vector dot product instead of the Euclidean distance ratio of the traditional method for matching. Perform a dot product operation on a certain feature point in the reference image and all the feature points in the image to be registered, then perform an inverse cosine transform on the dot product result and store it in the specified array, and then find the minimum angle θ1 and times Small angle θ2, if the ratio of the two is less than the specified threshold M, the feature point corresponding to the minimum angle is considered to match the feature point in the reference image. The reference image is a schematic diagram 3a of the feature point extraction result, and the image to be registered is another schematic diagram 3b of the feature point extraction result. As an example, M is 0.6. As shown in formula (7).

However, the traditional method is obtained by calculating the Euclidean distance between two feature point descriptors. Euclidean distance calculation formula is shown in formula (8).

In formula (8), Des _p and Des _q represent the feature descriptors of feature points p and q, respectively, and d is the Euclidean distance between the two feature points of p and q. For the feature descriptor e _l in the reference image, find the feature descriptors e _r and e _q closest to e _l in the image to be registered. Calculate the ratio N of the Euclidean distance D(e _l , e _r ) and D(e _l , e _q ), as shown in formula (9).

In formula (9), e _l = (e _l1 ,e _l2 ,…,e _l128 ), e _r =(e _r1 ,e _r2 ,…,e _r128 ), e _q =(e _q1 ,e _q2 ,…,e _q128 ). A threshold M is set, and the range of the threshold M is 0.4-0.6. If N<M, the pair of feature points (e _l , e _r ) is retained as a matching pair, otherwise discarded. Although the traditional method can find a suitable matching pair, when using the Euclidean distance ratio for calculation, it needs to perform square and square root operations many times, and the calculation complexity is high, which will increase the matching time and affect the matching efficiency. However, the calculation method for feature matching using the ratio of unit vector dot product arccosine angles adopted in the embodiment of the present invention only needs to perform basic operations such as vector multiplication and inverse cosine function, which simplifies the calculation formula and improves the matching efficiency of feature points .

As shown in FIG. 4 , it is a schematic diagram of a feature point matching result according to an embodiment of the present invention.

Step 5, removing mismatching points. In order to improve image registration accuracy, the embodiment of the present invention uses RANSAC algorithm to further purify feature points. Although the step 104 achieves a rough matching of feature points of the reference image and the image with registration, there are still mismatches. Therefore, the embodiment of the present invention uses the RANSAC algorithm to purify matching points. The specific process is as follows: randomly select 4 pairs from the obtained matching points, linearly calculate the transformation matrix H, calculate the distance d from each matching point to the corresponding matching point after H transformation, determine the distance threshold p, and satisfy d<p as the interior point, and re-estimate the transformation matrix according to the interior point. After N times of iterations, the initial value H _a of the transformation matrix H can be obtained after removing the mismatching points. The transformation matrix H is an 8-parameter 3×3 matrix, and project the transformation matrix H as shown in formula (10).

In formula (10), ( _xi , y _i ), (x' _i , y' _i ) are the corresponding matching points A _l and A _r of the reference image and the image to be stitched, respectively.

As shown in FIG. 5 , it is a schematic diagram of a matching result after removing mismatch points according to an embodiment of the present invention.

Step 6: Refining the transformation matrix. In the embodiment of the present invention, the L-M nonlinear iterative algorithm is used to refine the transformation matrix H, and the value of each parameter of the transformation matrix is continuously corrected through iterative calculation. In order to obtain the final accurate transformation matrix.

The flow chart of refining the transformation matrix by the L-M nonlinear iterative method is shown in Figure 6, and Figure 6 is a flow chart of refining the transformation matrix by the L-M nonlinear iterative method, including steps:

S601, set n=0; μ=0.01; set a threshold τ; the maximum number of iterations is N; the initial transformation matrix is m ₀ .

S602. Calculate the value of F(m _n ).

S603, judging whether F(m _n ) is greater than τ, if so, go to S604, otherwise, end the process.

S604, judge n! Whether it is greater than 0, if so, go to S605, otherwise, go to S606.

S605, judge whether F(m _n ) is smaller than F(m _n-1 ), if so, go to S606, otherwise, set m _n =m _n-1 , μ=10×μ and go to S607.

S606, set μ=μ/10 and calculate the value of J(m _n ), go to S607.

S607, judge whether n is greater than N, if so, end the process, otherwise, according to the formula:

m _n+1 ＝m _n －[J(m _n ) ^τ +J(m _n )μI] ^－1 J(m _n ) ^τ f(m _n ), where J(m) is the Jacobian matrix, let n =n+1, go to S602, that is, calculate the value of F(m _n ).

After the above steps, use the value of F(m _n ) into formula (11) - formula (13)

After calculation, the final transformation moment H _b is obtained. In the formula (m ₁₁ ~m ₃₂ ) are the components of the transformation matrix H, and in the formula (13), F(m) is the sum of the distances between all corresponding matching feature points. The goal of the LM algorithm is to solve the optimal solution of the eight parameters _m11 ~ _m32 of the transformation matrix, so that the sum of the distances between all corresponding matching feature points in the schematic diagram of the matching result after removing the mismatching points is the smallest, that is, When F(m) takes the minimum value, the corresponding solutions of the respective components m ₁₁ -m ₃₂ . According to the optimal solution of each component, the final transformation moment H _b can be obtained.

The image to be fused can be obtained after image transformation is performed on the original image of the image to be registered by using the final transformation matrix _Hb .

As shown in Fig. 7a, it is the original image of the image to be registered according to the embodiment of the present invention. After using the L-M nonlinear iterative algorithm to refine the transformation matrix from the original image of the image to be registered in the embodiment of the present invention in FIG. 7a, the image to be fused in the embodiment of the present invention is obtained, that is, FIG. 7b.

Step 7, image fusion, the embodiment of the present invention uses a weighted average fusion algorithm for image fusion. The weighted fusion algorithm is to directly take the same weight for the pixel values of the original image of the reference image and the corresponding pixel points of the image to be fused, and then perform a weighted average to obtain the pixel values of the corresponding pixel points of the fused image. As shown in FIG. 8 , it is the final spliced image according to the embodiment of the present invention.

As shown in Figure 9, it is a structural diagram of an aerial image mosaic system of an unmanned aerial vehicle forest area according to an embodiment of the present invention, including:

The image input module 801 is used to acquire aerial images of the forest area above the target area.

In the image enhancement module 802, due to the single color of the aerial image in some areas and no obvious change, the contrast of the aerial image is too low after being converted into a grayscale image, and the image details are blurred, so it is difficult to effectively use the SIFT algorithm to extract feature points. Therefore, the image enhancement module 802 is used to preprocess the image using the global histogram equalization algorithm, which can improve image contrast, increase image detail information, and thus increase the number of extracted feature points.

The specific processing process is as follows: normalize the gray level r of the image _gzh(x,y) to the interval [0,1], when r=0 it is black, and when r=1 it is white. g _{zh(x, y)} gray scale range is [0, L-1], the total number of pixels is n, then the number of pixels with gray scale r _k is n _k . The transformation corresponding to the global histogram equalization is shown in formula (1).

In formula (1), P _r (r _j ) is the probability distribution function of g _zh(x,y) , and S _k is the gray level k of the image g _zh(x,y) from 0 to (L-1) When , the probability distribution function P _r (r _j ) of the image g _zh(x,y) is summed, L is the total number of gray levels of the image, and T(r) represents the transformation function. Using the transformation formula to change the grayscale of the aerial image, the contrast effect after image enhancement is obtained as shown in Figures 2a-2c, the contrast of the enhanced grayscale aerial image is obviously improved, and the image detail information is more abundant. 2a is the original aerial image of the embodiment of the present invention, FIG. 2b is the grayscale aerial image of the embodiment of the present invention, and FIG. 2c is the enhanced grayscale aerial image of the embodiment of the present invention.

Feature point extraction module 803. When the UAV is shooting the image of the target area, due to the interference of illumination changes and other environmental factors, the input image sequence may be rotated, and the viewing angle and scale of the image may also change. In order to overcome the above interference , the embodiment of the present invention uses the SIFT algorithm to perform image registration. The image registration is the core step of image stitching, and the accuracy of registration will directly affect the final effect of image stitching. The image registration is to calculate the space transformation parameters between the reference image and the image to be registered, determine the space transformation model, and realize the transformation of the coordinate system between the two images. The feature point extraction module 803 is used to extract feature points in scale space using the SIFT algorithm, which has good scale invariance and robustness to changes in viewing angles.

The specific functions of the feature point extraction module 803 include:

Scale space extreme value detection, the scale space L(x,y,σ) of a two-dimensional image is defined as the convolution of the scale-variable Gaussian function G(x,y,σ) and the original image I(x,y), As shown in formulas (2) and (3).

L(x,y,σ)=G(x,y,σ)*I(x,y) (3)

In formula (3), * is the convolution operation, (x, y) is the position of the image pixel, σ is the scale space factor, the larger the value of σ, the blurrier the image, and the smaller the value of σ, the clearer the image. In order to effectively detect relatively stable feature points in the scale space, the Gaussian difference scale space D(x, y, σ) is generated by convolution of Gaussian difference kernels of different scales with the image. As shown in formula (4).

D(x,y,σ)=(G(x,y,kσ)-G(x,y,σ))*I(x,y)=L(x,y,kσ)-L(x,y ,σ)

(4)

In formula (4), k is a constant multiplicative factor, and its value is related to the number S of pictures per layer in the scale space, that is, k=2 ^1/S . The enhanced grayscale aerial image is detected, and if the grayscale value of a certain pixel reaches an extreme value in its surrounding and 26 neighborhoods, record the position and scale of this point as a candidate feature point. This patent combines the contrast threshold screening method with the features and characteristics of forest images to detect extreme points, reduce the extraction of invalid feature points in the image, and improve the efficiency of the algorithm for detecting and extracting feature points of this type of image. The principle of threshold selection is to ensure the correct registration in the forest image while minimizing feature point extraction. The relationship between the threshold size and the correct registration rate of feature points is obtained through experimental verification, so as to select the optimal contrast threshold. The main basis of this method is: there are a large number of similar structures in the forest image, and the pixel value changes between the forest scenes are small, but there are obvious contrast and pixel value changes between the forest and non-forest in the image. Therefore, the method of increasing the contrast threshold can be used to extract the main scene in the image content, thereby extracting feature points that are more likely to be useful for image registration, and reducing the detection of feature points that do not contribute to registration.

Precisely locate feature points, accurately determine the position and scale of feature points by fitting three-dimensional quadratic functions, and filter low-contrast feature points and unstable edge response points to enhance matching stability and improve noise resistance.

The main direction of the feature point is determined. In order to make the descriptor invariant to rotation, a direction is assigned to each feature point by using the local features of the grayscale aerial image. Using the characteristics of the gradient and direction distribution of the neighborhood pixels of the feature point, the modulus m(x, y) and direction (x, y) of the gradient are obtained, as shown in formulas (5) and (6).

To generate a feature point descriptor, in order to ensure the invariance of feature vector rotation, first rotate the coordinate axis to the main direction of the feature point, use a Gaussian circular window to Gaussian weight the gradient modulus, and use the feature point as the center of the 16×16 window, Divide this window into 4×4 sub-regions, each region counts 8 directions through the gradient histogram, and finally forms a 128-dimensional SIFT feature vector that is not affected by scale, rotation and illumination.

Feature points are determined using the feature point descriptor.

As shown in Figure 3a and Figure 3b, it is the feature point extraction result. Fig. 3a is a schematic diagram of a feature point extraction result of an embodiment of the present invention, and Fig. 3b is a schematic diagram of a feature point extraction result of another embodiment of the present invention.

The feature point matching module 804 is used to match the feature points. In order to simplify the matching process and improve the feature point matching speed, the feature point matching module 804 uses the ratio of the arccosine function of the unit vector dot product to replace the Euclidean distance of the traditional method The ratio is matched. Perform a dot product operation on a certain feature point in the reference image and all the feature points in the image to be registered, then perform an inverse cosine transform on the dot product result and store it in the specified array, and then find the minimum angle θ1 and times Small angle θ2, if the ratio of the two is less than the specified threshold M, the feature point corresponding to the minimum angle is considered to match the feature point in the reference image. The reference image is a schematic diagram 3a of the feature point extraction result, and the image to be registered is another schematic diagram 3b of the feature point extraction result. As an example, M is 0.6. As shown in formula (7).

However, the traditional method is obtained by calculating the Euclidean distance between two feature point descriptors. Euclidean distance calculation formula is shown in formula (8).

In formula (8), Des _p and Des _q represent the feature descriptors of feature points p and q, respectively, and d is the Euclidean distance between the two feature points of p and q. For the feature descriptor e _l in the reference image, find the feature descriptors e _r and e _q closest to e _l in the image to be registered. Calculate the ratio N of the Euclidean distance D(e _l , e _r ) and D(e _l , e _q ), as shown in formula (9).

In formula (9), e _l = (e _l1 ,e _l2 ,…,e _l128 ), e _r =(e _r1 ,e _r2 ,…,e _r128 ), e _q =(e _q1 ,e _q2 ,…,e _q128 ). A threshold M is set, and the range of the threshold M is 0.4-0.6. If N<M, the pair of feature points (e _l , e _r ) is retained as a matching pair, otherwise discarded. Although the traditional method can find a suitable matching pair, when using the Euclidean distance ratio for calculation, it needs to perform square and square root operations many times, and the calculation complexity is high, which will increase the matching time and affect the matching efficiency. The feature point matching module 804 uses the ratio of the unit vector dot product arccosine angle to perform feature matching. The calculation method only needs to perform basic operations such as vector multiplication and inverse cosine function, which simplifies the calculation formula and improves the matching efficiency of feature points. .

As shown in FIG. 4 , it is a schematic diagram of a feature point matching result according to an embodiment of the present invention.

The module 805 for removing incorrectly matched points is used to remove incorrectly matched feature points. In order to improve the accuracy of image registration, the embodiment of the present invention uses the RANSAC algorithm to further purify the feature points. Although the feature point matching module 804 realizes the rough matching of the feature points of the front and back images, there are still mismatches. Therefore, the module 805 for removing incorrect matching points uses a RANSAC (Random Sample Consensus, Random Sample Consensus) algorithm to purify matching points. The specific process is as follows: randomly select 4 pairs from the obtained matching points, linearly calculate the transformation matrix H, calculate the distance d from each matching point to the corresponding matching point after H transformation, determine the distance threshold p, and satisfy d<p as the interior point, and re-estimate the transformation matrix according to the interior point. After N iterations, the initial value Ha of the transformation matrix H can be obtained. The transformation matrix H is a 3×3 matrix containing 8 parameters, and the projected The above transformation matrix H is shown in formula (10).

In formula (10), (x _i , y _i ), (x′ _i , y′ _i ) are the corresponding matching points A _l and A _r of the reference image and the image to be stitched, respectively.

As shown in FIG. 5 , it is a schematic diagram of a matching result after removing mismatch points according to an embodiment of the present invention.

The refining transformation matrix module 806 is configured to refine the transformation matrix H by using the L-M nonlinear iterative algorithm.

The flow chart of refining the transformation matrix by the L-M nonlinear iterative method is shown in Figure 6, and Figure 6 is a flow chart of refining the transformation matrix by the L-M nonlinear iterative method, including steps:

S601, set n=0; μ=0.01; set a threshold τ; the maximum number of iterations is N; the initial transformation matrix is m ₀ .

S602. Calculate the value of F(m _n ).

S603, judging whether F(m _n ) is greater than τ, if so, go to S604, otherwise, end the process.

S604, judge n! Whether it is greater than 0, if so, go to S605, otherwise, go to S606.

S605, judge whether F(m _n ) is smaller than F(m _n-1 ), if so, go to S606, otherwise, set m _n =m _n-1 , μ=10×μ and go to S607.

S606, set μ=μ/10 and calculate the value of J(m _n ), go to S607.

S607, judge whether n is greater than N, if so, end the process, otherwise, according to the formula:

m _n+1 ＝m _n －[J(m _n )τ+J(m _n )μI] ^－1 J(m _n ) ^τ f(m _n ), where J(m) is the Jacobian matrix, let n =n+1, go to S602, that is, calculate the value of F(m _n ).

After the above steps, use the value of F(m _n ) into formula (11) - formula (13)

After calculation, the final transformation moment H _b is obtained. In the formula (m ₁₁ ~m ₃₂ ) are the components of the transformation matrix H, and in the formula (13), F(m) is the sum of the distances between all corresponding matching feature points. The goal of the LM algorithm is to solve the optimal solution of the eight parameters _m11 ~ _m32 of the transformation matrix, so that the sum of the distances between all corresponding matching feature points in the schematic diagram of the matching result after removing the mismatching points is the smallest, that is, When F(m) takes the minimum value, the corresponding solutions of the respective components m ₁₁ -m ₃₂ . According to the optimal solution of each component, the final transformation moment H _b can be obtained.

The image to be fused can be obtained by performing image transformation on the image to be registered using the final transformation matrix _Hb .

As shown in Fig. 7a, it is the original image of the image to be registered according to the embodiment of the present invention. After using the L-M nonlinear iterative algorithm to refine the transformation matrix from the original image of the image to be registered in the embodiment of the present invention in FIG. 7a, the image to be fused in the embodiment of the present invention is obtained, that is, FIG. 7b.

The image fusion module 807 is configured to perform image fusion using a weighted average fusion algorithm. The weighted fusion algorithm is to directly take the same weight for the pixel values of the original image of the reference image and the corresponding pixel points of the image to be fused, and then perform a weighted average to obtain the pixel values of the corresponding pixel points of the fused image. As shown in FIG. 8 , it is the final fused image of the embodiment of the present invention.

To sum up, the present invention provides a method and system for splicing aerial images of unmanned aerial vehicles in forest areas. First, the global histogram equalization algorithm is used to preprocess the images, and then the SIFT feature image splicing algorithm is used to extract feature points. After removing the mismatched feature points, the transformation matrix is refined, and finally the image is fused to obtain the final fused image. The method and system for splicing aerial images of unmanned aerial vehicle forest areas have the characteristics of fast speed and high accuracy, and can quickly splice aerial images with single color and low contrast, which is convenient for users to analyze subsequent remote sensing images and processing.

Those of ordinary skill in the art should understand that: the discussion of any of the above embodiments is exemplary only, and is not intended to imply that the scope of the present disclosure (including claims) is limited to these examples; under the idea of the present invention, the above embodiments or Combinations between technical features in different embodiments are also possible, steps may be carried out in any order, and there are many other variations of the different aspects of the invention as described above, which are not presented in detail for the sake of brevity. Therefore, any omissions, modifications, equivalent replacements, improvements, etc. within the spirit and principles of the present invention shall be included within the protection scope of the present invention.

Claims (7)

1. A kind of unmanned aerial vehicle forest area aerial photography image splicing method is characterized in that, comprises the following steps: Take images of the forest area over the target area through the drone equipped with visible light acquisition equipment; Using a global histogram equalization algorithm to preprocess the aerial image to obtain an enhanced grayscale aerial image; Feature point detection and extraction are performed on the enhanced grayscale aerial image through an optimized contrast threshold screening method, and an improved feature point matching method is used to match the feature points of the reference image and the image to be registered to obtain the feature Point matching diagram; Use the RANSAC algorithm to remove mismatching points and calculate the initial value of the transformation matrix, then use the L-M nonlinear iterative algorithm to refine the calculated initial transformation matrix, and continuously correct the values of the parameters of the transformation matrix through iterative calculations to obtain the final accurate value Transformation matrix, using the final transformation matrix to obtain the image to be fused after performing image transformation on the original image of the image to be registered; performing image fusion on the original image of the reference image and the image to be fused by using a weighted average fusion algorithm to obtain a final fused image. 2. a kind of unmanned aerial vehicle forest region aerial photography image splicing method according to claim 1, is characterized in that, by unmanned aerial vehicle according to the mode that plan route scans back and forth to photograph image, the overlapping degree between the image of described photographed image is 40%-60%. 3. a kind of unmanned aerial vehicle forest region aerial image mosaic method according to claim 1, is characterized in that, described global histogram equalization algorithm processing process is as follows: the gray level of aerial image g _{zh (x, y)} Level r is normalized to the interval [0,1], black when r=0, white when r=1; _gzh(x,y) grayscale range is [0,L-1], the total number of pixels is n, then the number of pixels with gray level r _k is n _k ; the transformation corresponding to the global histogram equalization is shown in formula (1); P _r (r _j ) is the probability distribution function of g _{zh(x, y)} , and S _k is when the gray level k of the aerial image g _{zh(x, y)} is taken from 0 to (L-1), Summing the image _{gzh(x, y)} probability distribution function P _r (r _j ), L is the total number of gray levels of the aerial image, T (r) represents a transformation function; use the transformation formula to change the aerial image The grayscale of the enhanced grayscale aerial image is obtained. 4. a kind of unmanned aerial vehicle forest region aerial image splicing method according to claim 1, it is characterized in that, the contrast threshold screening method of described optimization is: first remove the point that pixel value in image is lower than contrast threshold, described The threshold selection method is as follows: the relationship between the threshold size and the correct registration rate of feature points is obtained through experimental verification, and the optimal contrast threshold is selected. 5. a kind of unmanned aerial vehicle forest area aerial photography image splicing method according to claim 1, is characterized in that, described improved feature point matching algorithm processing procedure is as follows: a certain feature point in the reference image and the image to be spliced Do the dot product operation for all the feature points in , and then store the dot product result in the specified array by inverse cosine transform, and then find the minimum angle θ ₁ and the second smallest angle θ ₂ from the array, if the ratio of the two is less than the specified Threshold M=0.6, then it is considered that the feature point e _l corresponding to the minimum angle matches the feature point e _r in the reference image; (e _l , e _r ) is left as a matching pair of feature points, and the matching formula is as in formula (2) Shown: In this way, initial matching can be performed on the extracted feature points to obtain initial matching point pairs. 6. a kind of unmanned aerial vehicle forest region aerial photography image splicing method according to claim 1, it is characterized in that, RANSAC algorithm and LM non-linear iterative algorithm are combined, and described RANSAC algorithm removes the concrete processing process of mismatch point as follows: Randomly select 4 pairs from the obtained matching points, linearly calculate the transformation matrix H, calculate the distance d from each matching point to the corresponding matching point after H transformation, determine the distance threshold p, and use the point satisfying d<p as the interior point , and re-estimate the transformation matrix according to the internal points. After N iterations, the initial value H _a of the transformation matrix H can be obtained after removing the mismatch points. The formula for removing the mismatch points by the RANSAC algorithm is shown in formula (3): In formula (3), H is a transformation matrix, and the transformation matrix H is a 3×3 matrix containing 8 unknown parameters (m ₁₁ ～m ₃₂ ); ( _xi , y _i ), (x′ _i , y ' _i ) are corresponding matching points A _l and A _r in the reference image and the image to be stitched, respectively. The initial value of the calculated transformation matrix is H _a . The method of refining the calculated initial transformation matrix H _a by using the LM nonlinear iterative algorithm, continuously correcting the values of the parameters of the transformation matrix through iterative calculations, and obtaining the final accurate transformation matrix H _b includes the steps of: S601, set n=0; μ=0.01; set a threshold τ; the maximum number of iterations is N; the initial transformation matrix is m ₀ ; S602, calculating the value of F(m _n ); S603, judge whether F(m _n ) is greater than τ, if so, go to S604, otherwise end the process; S604, judge n! Whether it is greater than 0, if so, go to S605, otherwise go to S606; S605, judge whether F(m _n ) is smaller than F(m _n-1 ), if so, go to S606, otherwise, set m _n =m _n-1 , μ=10×μ and go to S607; S606, set μ=μ/10 and calculate the value of J(m _n ), go to S607; S607, judge whether n is greater than N, if so, end the process, otherwise, according to the formula: m _n+1 ＝m _n －[J(m _n ) ^τ +J(m _n )μI] ^－1 J(m _n ) ^τ f(m _n ), where J(m) is the Jacobian matrix, let n =n+1, go to S602, namely calculate the value of F(m _n ); Use the value of F(m _n ) into (11)-Equation (13) In the formula (m ₁₁ ~m ₃₂ ) are the components of the transformation matrix H, in the formula (13), F(m) is the sum of the distances between all corresponding matching feature points; to calculate the optimal solution of each component, you can The final transformation moment H _b is obtained. 7. A system for mosaicing aerial images of unmanned aerial vehicle forest areas, characterized in that it comprises: The image input module is used to obtain the aerial photograph image over the target forest area; An image enhancement module, configured to use a global histogram equalization algorithm to preprocess the aerial image to obtain an enhanced grayscale aerial image; The feature point extraction module is used to use the optimized contrast threshold screening method to perform scale space extremum detection on the enhanced gray-scale aerial image, accurately locate feature points, determine the main direction of feature points, and generate feature point descriptors; use the The feature point descriptor is used to determine the feature point, and the reference image and the image to be registered of the feature point are obtained; The feature point matching module is used to perform feature point matching on the feature points of the reference image and the image to be registered using the ratio of the arccosine function of the unit vector dot product to obtain a schematic diagram of feature point matching; Remove the mismatching point module, which is used to make the RANSAC algorithm remove the mismatching points in the feature point matching schematic diagram, obtain the matching result schematic diagram after removing the mismatching points and calculate and obtain the initial transformation matrix; The refining transformation matrix module is used to refine the initial transformation matrix using the L-M nonlinear iterative algorithm, so that the sum of the distances between all corresponding matching feature points in the matching result schematic diagram after removing the mismatching points is the smallest, and obtains the final transformation matrix; Using the final transformation matrix to obtain the image to be fused after performing image transformation on the original image of the image to be registered; An image fusion module, configured to use a weighted average fusion algorithm to perform image fusion on the original image of the reference image and the image to be fused to obtain a final fused image.