CN106971404A - A kind of robust SURF unmanned planes Color Remote Sensing Image method for registering - Google Patents

Abstract

The invention relates to a robust SURF unmanned aerial vehicle color remote sensing image registration method, that is, in the feature point description stage, the feature point itself is transformed into a color space to obtain the color information of the feature point itself, and then expressed in the description operator . The matching works well and there are no cross lines appearing in the image.

Description

A Robust SURF UAV Color Remote Sensing Image Registration Method

technical field

The invention relates to a UAV image registration technology, in particular to a remote sensing image registration method.

Background technique

SURF (Speeded Up Robust Features) algorithm is a robust local feature detection algorithm. This algorithm not only has a highly repeatable detector and a highly distinguishable descriptor, but also has strong robustness and comparative High computing speed; the accuracy can reach the sub-pixel level, and it can ensure that the feature descriptor extracted from the image has good invariance under the conditions of scale scaling, brightness change, rotation, occlusion and noise. Currently, the algorithm It has been widely used in image registration technology. However, due to the small size and large number of UAV remote sensing images, the requirements for registration accuracy are extremely stringent. How to enhance the salience of color image feature points, so as to further improve the registration accuracy, especially the precise registration of color remote sensing images with high feature similarity has always been an urgent problem to be solved.

Some researchers have proposed the method of normalizing the intensity of RGB components and then using the SURF algorithm to process color remote sensing images and the method of color invariant features based on SURF, but the first method has limitations in describing color invariant features The color information cannot be fully utilized, and the second method takes a long time and has poor real-time performance, which cannot meet the actual requirements.

Contents of the invention

technical problem to be solved

Aiming at the problem that the registration accuracy is reduced by only using image brightness information and ignoring the color information when the color remote sensing image of the UAV is registered using the SURF algorithm, the present invention proposes a robust SURF UAV color remote sensing image registration method.

Technical solutions

A method for registration of robust SURF UAV color remote sensing images, characterized in that the steps are as follows:

Step 1: Perform grayscale processing on the acquired UAV color remote sensing reference image image1 and image2 to be registered to become a grayscale image, and perform Gaussian smoothing on the grayscale image to complete the scale space construction, and use the feature detection operator Hessian matrix The determinant of the feature points is extracted, and all the feature points of image1 and image2 are recorded as the target set {p _i }, i=1,2,...,n and the reference set {q _j },j=1,2, ...,m;

Step 2: For the obtained target set {p _i } and reference set {q _j }, a square area with a side length of 20s is determined centering on each feature point, where s represents the sampling interval in different scale spaces, and the surrounding Neighborhood information is divided into 4×4 sub-regions, and each small region is divided into 5×5 sampling points. Finally, Harr wavelet is used to calculate the response and modulus dx,|dx|,dy, |dy|, and counting the total response ∑dx, ∑|dx|, ∑dy, ∑|dy| of 5×5 sampling points, the 64-dimensional SURF feature descriptor of the feature point can be obtained:

V=(∑dx,∑|dx|,∑dy,∑|dy|)

Then the SURF descriptor of the feature point is defined in the horizontal or vertical direction as:

V _gray ＝(V ₁ ,V ₂ ,V ₃ ,V ₄ ,...,V ₁₆ )

Among them, V _i , _{i =} 1, 2, 3,..., 15, 16 represent the description vectors corresponding to the 16 sub-regions; from this we can get all the 64-dimensional SURF feature descriptor for feature points;

Step 3: For the feature point with coordinates (x, y), it is expressed as R _{(x, y)} , G _{(x, y)} , B _{(x, y)} in the RGB color space, and it is converted into the color space, and the conversion To the YIQ color space is Y _(x,y) , I _(x,y) , Q _(x,y) ;

The transformation from RGB color space to YIQ color space is expressed as:

Adding it to the original algorithm descriptor yields a new color descriptor as follows:

Normalize the above formula:

where _{i k} , _k =1,2,3,...,63,64 represent the 64-dimensional description vector of feature points in the classic SURF algorithm; 67-dimensional description vector of each feature point;

Step 4: Two-way matching of feature points

Step 4a: use the Euclidean distance formula to calculate the nearest Euclidean distance d _nt and the next closest Euclidean distance d _snt in the forward direction of the feature point p1 in the reference image image1;

Use the Euclidean distance formula to calculate the feature points of the feature point p1 in the reference image image1 in the forward direction of the nearest Euclidean distance d1 _nt and the second closest Euclidean distance d1 _snt in the image2, respectively, and the feature points are p2 _nt and p2 _snt ;

Step 4b: Calculate the ratio T ₁ of the closest distance d1 _nt and the second closest distance d1 _snt , namely

Compare T ₁ with the judgment threshold T: if T ₁ <T, proceed to step 4c, otherwise, jump out automatically;

Step 4c: For the feature point corresponding to p2 _nt , use the Euclidean distance formula to calculate the feature points of the reverse closest distance d2 _nt and the next closest distance d2 _snt in image1 are p1 _nt and p1 _snt respectively;

Step 4d: Calculate the ratio T ₂ of the closest distance to the next closest distance

Compare T ₁ with the judgment threshold T: if T ₂ <T, and the feature point p1 _nt is looking for is the same feature point as the original feature point p1, it is considered to be a correct matching point pair, otherwise return to step 4a for the reference image Other feature points in image1 are used to judge new feature points;

Traverse all the feature points in the reference image target set in turn, and finally achieve two-way matching and purification;

Step 5: Use the RANSAC algorithm to obtain the parameters of the transformation matrix

Step 5a: Find the feature point pair m _i '←→m _i after the two-way matching and purification after step 4 is completed, i=1,2,...n;

Step 5b: Randomly select 4 pairs of matching points from the set of matching points, and use the coordinate values of these 4 pairs of matching points to obtain the transformation matrix H;

Step 5c: Based on the Euclidean distance, find the point pair in the matching point pair set that meets the condition of d(H _i m _i ,m _i ')<t, where H _i m _i represents the matrix of the feature points in image2 The transformation is mapped to the position coordinates in image1, m _i 'indicates the position coordinates of _mi in image1 corresponding to the feature point, d(H _i m _i ,m _i ') indicates the Euclidean distance between two coordinate vectors, and will meet the conditions The point pair is used as the final interior point, and the number of interior points satisfying the H _i constraint is recorded;

Step 5d: Repeat steps 5b and 5c n times, and record the number of interior points each time;

Step 5e: Select H _best corresponding to the largest number of interior points, find all matching point pairs that meet the condition of d(H _best m _i ,m _i ')<t, and use them as the final interior points, that is, the correct matching point pairs, Wrong matching point pairs that do not meet the condition of d(H _best m _i ,m _i ')<t are outliers and should be eliminated;

Step 5f: According to the random sampling consensus algorithm, N matching point pairs are obtained, and the 2N matrices composed of x ₁ , y ₁ , x ₂ , y ₂ are marked, denoted as A, and the SVD singular value is performed on it Decompose A=UDV ^T , U and V are the singular vectors of A, D is the singular values on the diagonal in descending order, and the last column of V is reconstructed into a 3*3 matrix, which is the required perspective transformation matrix Parameter Estimation;

Step 6: Bicubic interpolation for registration

If the coordinates of any point on the reference image image1 are (x, y), its gray value is g(x, y), and the corresponding point coordinates of this point on the image to be registered are (u, v), and its gray value The value is f(u,v), and [ ] means taking an integer, and the matrix composed of 16 pixels in the 4*4 neighborhood of the corresponding point in the image to be registered is B:

The function f(u+i,v+j) in the elements in the matrix B, i=-1,0,1,2; j=-1,0,1,2, defined as the corresponding point on the image to be registered The gray value when the coordinates are (u+i,v+j);

Then the calculation formula of the gray value of the point to be interpolated in the reference image is as follows:

g(x,y)=f(u,v)=ABC

in:

s(w) is the basis function of bicubic interpolation, which is an approximation to sin(x*π)/x. Through bicubic interpolation, the image2 image is interpolated into the reference image image1 to achieve the final high-precision registration.

The judgment threshold T in the step 4b is 0.4-0.8.

t in the step 5c is selected as 4.

100≤n≤200 in the step 5d.

Beneficial effect

A robust SURF UAV color remote sensing image registration method proposed by the present invention, that is, in the feature point description stage, the color space transformation is performed on the feature point itself to obtain the color information of the feature point itself, and then in the description operator. expression. The matching works well and there are no cross lines appearing in the image.

Description of drawings

Fig. 1 is the color remote sensing images image1 and image2 of the UAV in the present invention, wherein Fig. 1(a) is the reference image image1, and Fig. 1(b) is the image to be registered image2.

Figure 2 shows the registration results of the original SURF algorithm.

Fig. 3 is a robust SURF registration result of the present invention.

detailed description

Now in conjunction with embodiment, accompanying drawing, the present invention will be further described:

Because the components of the RGB color space are closely related to the brightness, and the color remote sensing images usually collected are too sensitive to the brightness, the RGB color space is not suitable for color image processing. The YIQ color space (the Y component represents the brightness information of the image, and the I and Q two components carry the color information) and the RGB color space can be directly transformed through a matrix multiplication, and it is a linear transformation with high real-time performance. On the other hand, YIQ has the ability to extract the brightness component, and the color information is also separated at the same time, with a small amount of calculation and good clustering characteristics, so it can be effectively used for color image processing. This method transforms the feature points from the RGB space to the YIQ space in the feature point description stage, expresses the feature points in a YIQ model, adds the color information of the feature points to the descriptor, and adds Three-dimensional color information vectors constitute 67-dimensional descriptors to enhance the discrimination of color image feature point descriptors. Then, RANSAC (Random Sample Consensus), bidirectional matching, and bicubic interpolation algorithms are combined to achieve precise registration of UAV color remote sensing images. It mainly includes the following steps:

Step 1. Use the SURF algorithm to perform feature point detection to obtain the target set and reference set

Let the UAV color remote sensing image image1 be the reference image, the feature point description subset in it is the target set {p _i }, i=1,2,...,n, another color remote sensing image image2 is the image to be registered , where the feature point description subset is the reference set {q _j }, j=1,2,...,m;

Gray-scale processing is performed on the acquired UAV color remote sensing images image1 and image2 to become gray-scale images, Gaussian smoothing processing is performed to complete the construction of the scale space, and feature point extraction is performed using the determinant of the feature detection operator Hessian matrix to obtain image1, All feature points of image2 are respectively recorded as the target set {p _i }, i=1,2,...,n and the reference set {q _j },j=1,2,...,m, the specific description is as follows:

That is, for a certain point X=(x,y) in the grayscale image I(x,y) obtained after grayscale processing, the Hessian matrix on the scale σ is defined as:

where: L _xx (X,σ) is the Gaussian second order differential with scale σ Convolved with image I(x,y) at point X:

in, Represents a two-dimensional convolution operation. Similarly, L _xy (X, σ) and L _yy (X, σ) are Gaussian second-order differentials with a scale of σ Convolved with I(x,y) at point X.

Determinant of the Hessian matrix:

Det(H)＝L _xx (X,σ)L _yy (X,σ)-L ² _xy (X,σ) (3)

Different scales σ represent different spatial resolutions, different scales of σ values have different images after Gaussian smoothing, and the scale space is constructed:

Between-group scale:

σ _octave = σ ₀ *2 ^octave-1 (4)

Intragroup scale:

σ-----scale space coordinate, which is a numerical representation of scale space;

s ----- scale space coordinates within the group, which is a numerical representation of the scale space within the group;

σ ₀ -----Initial scale, determined according to camera performance, generally choose 1.6;

S-----the number of layers per group (4 layers);

octave-----number of groups (4 groups);

The scale space is established by formulas (4) and (5), and then the local extremum of the determinant response is searched in the scale space of the image through the three-dimensional non-maximum suppression search technology to complete the detection of feature points;

Step 2. Description of feature points and adding color information

2.1 Feature point description

For the obtained target set {p _i } and reference set {q _j }, a square area with a side length of 20s (s represents the sampling interval in different scale spaces) is determined with each feature point as the center, and the surrounding neighborhood information Divide into 4×4 sub-regions, each small region is divided into 5×5 sampling points, and finally use Harr wavelet to calculate the vertical and horizontal response and modulus dx,|dx|,dy,|dy| of each small region , and counting the total response ∑dx, ∑|dx|, ∑dy, ∑|dy| of 5×5 sampling points, the 64 (4×4×4) dimensional SURF feature descriptor of the feature point can be obtained:

V=(∑dx,∑|dx|,∑dy,∑|dy|) (6)

Assuming that the coordinates of a feature point are (x, y), the SURF descriptor of the feature point can be defined in one direction (horizontal or vertical) as:

V _gray ＝(V ₁ ,V ₂ ,V ₃ ,V ₄ ,...,V ₁₆ ) (7)

Wherein, V _i (i=1, 2, 3, . . . , 15, 16) represents descriptors corresponding to 16 sub-regions. The 64-dimensional SURF feature descriptors of all the feature points in the target set {p _i } and reference set {q _j } can be obtained by using the above method.

2.2 Add color information

The color information of the feature point itself is added to the original descriptor, and the coordinates are (x, y) feature points, which are expressed as R _{(x, y)} , G _{(x, y)} , B _(x , y) in the RGB color space _y) , convert it to color space, and convert it to YIQ color space as Y _(x,y) , I _(x,y) , Q _(x,y) , RGB color space is most commonly used in daily life. However, the color remote sensing images we usually collect are too sensitive to brightness, and the components of RGB color space are closely related to the brightness. Therefore, the RGB color space is not suitable for color image processing. However, YIQ and RGB can be directly transformed into each other through a matrix multiplication, and become a linear transformation with high real-time performance. On the other hand, the YIQ color space can extract the luminance component and separate the color information at the same time, the calculation amount is small, and the clustering characteristics are good, so it can be effectively used for color image processing, which is the root of our improvement.

The transformation from RGB color space to YIQ color space is expressed as:

Adding it to the original algorithm descriptor yields a new color descriptor as follows:

V _color ＝(V ₁ ,V ₂ ,V ₃ ,V ₄ ,...,V ₁₆ ,Y _(x,y) ,I _(x,y) ,Q _(x,y) ) (9)

In order to make the expression clearer, the simple transformation of formula (6)

V=(i ₁ ,i ₂ ,i ₃ ,i ₄ ) (10)

Then formula (9) is redefined as follows:

V _color ＝(i ₁ ,i ₂ ,i ₃ ,i ₄ ,...,i ₆₄ ,Y _(x,y) ,I _(x,y) ,Q _(x,y) ) (11)

For better robustness to rotation, scale, lighting. Here, formula (11) needs to be normalized:

Where |V _color | is the calculation mode of formula (11), the expression is:

Among them, i _k (k=1,2,3,...,63,64) represents the 64-dimensional descriptor of the feature point in the classical SURF algorithm. From formula (11), we can see that on the basis of the original 64-dimensional feature point descriptor, the color information of the feature point itself is also loaded in to form a 67-dimensional descriptor. Compared with the original algorithm , there is not much change in time, but in terms of performance, the accuracy will be further improved due to the addition of color information.

By analogy, the 67-dimensional descriptors of each feature point in the target set {p _i } and reference set {q _j } are finally obtained.

Step 3. Feature point matching

For each feature point in the target set in the reference image image1, query its nearest neighbor Euclidean distance d _nt and second nearest neighbor Euclidean distance d _snt in the reference set, if it satisfies:

Then retain the matching between the feature points in the target set conforming to the formula (14) and the nearest neighbors of the feature points in the target set that are queried in the benchmark set; otherwise, this matching pair is eliminated. T in the formula is the judgment threshold, Take a value between 0.4-0.8.

Through the above steps, we get the matching point pairs under the ratio measurement criterion of the nearest neighbor Euclidean distance and the second nearest neighbor Euclidean distance. These matching point pairs are still highly likely to be mis-matched, and need to be purified. Here, the two-way matching method is used, and the operation is as follows:

To make the expression clearer, the descriptor formula (12) in image1 and image2 is recorded here as v1=(v1 ₁ ,v1 ₂ ,v1 ₃ ,v1 ₄ ,...,v1 ₆₇ ) and v2=(v2 ₁ , v2 ₂ ,v2 ₃ ,v2 ₄ ,...,v2 ₆₇ ), from the Euclidean distance formula can be obtained:

Among them, the Euclidean distance between the two vectors v1 and v2 is expressed as d(v1,v2), v1 _i , v2 _i (i=1,2,3,...,66,67) represent the vectors v1, v2 Each component value of ;

The closest distance and the second closest distance are denoted as d _nt and d _snt , and T _i is the ratio of the closest distance to the second closest distance, expressed as follows:

According to the formulas (15) and (16), the two-way matching point pair is purified, and the specific steps are as follows:

(1) Referring to the feature point p1 in image1, use the Euclidean distance formula to find the feature points of the closest Euclidean distance and the next closest Euclidean distance in image2, which are p2 _nt and p2 _snt respectively, and use d1 _nt and d1 _snt to describe the positive The shortest distance and the second shortest distance in the direction, where the direction from the feature point of the reference image target set in the reference set is located, and the forward direction is determined, that is, a feature point in the target set is traversed in the reference set to find a point that satisfies the condition; and from the reference set The direction to look for in the target set is defined as reverse;

(2) Use formula (16) to calculate the ratio T ₁ of the closest distance to the second closest distance, namely

If T ₁ <T, proceed to step (3), otherwise, jump out automatically;

(3) For the feature point corresponding to p2 _nt , use the formula (15) to find the feature points with the shortest distance and the second shortest distance in image1, which are p1 _nt and p1 _snt respectively, and use d2 _nt and d2 _snt to describe the reverse shortest distance and the next closest distance;

(4) Use formula (16) to calculate the ratio T ₂ of the closest distance and the next closest distance

If T ₂ <T, and at the same time the feature point that p1 _nt is looking for is the same feature point as the original feature point p1, then it is considered to be a correct matching point pair, otherwise return to step (1) to update other feature points in the reference image image1 feature point judgment;

It traverses all the feature points in the target set of the reference image in turn, and finally achieves two-way matching and purification.

Step 4, RANSAC algorithm to obtain transformation matrix parameters

Firstly, solve the transformation matrix. Given two points m' and m on the color remote sensing images image1 and image2, according to the perspective transformation matrix H, there are: m'～Hm, where ~ means equal in proportion, and the coordinates of m and m' are respectively (x ₁ ,y ₁ ) and (x ₂ ,y ₂ ), written in homogeneous coordinate form: (x ₁ ,y ₁ ,z ₁ ) and (x ₂ ',y ₂ ',z ₂ ), where x ₂ '/z ₂ =x ₂ , y ₂ '/z ₂ =y ₂ , then the formula

Then get (20) and (21) on the basis of homogeneous coordinates:

Without loss of generality and for the convenience of calculation, let z ₁ =1, from (20) and (21):

x ₁ H ₁₁ +y ₁ H ₁₂ +H ₁₃ -x ₁ x ₂ H ₃₁ -y ₁ x ₂ H ₃₂ -x ₂ H ₃₃ =0 (22)

x ₁ H ₂₁ +y ₁ H ₂₂ +H ₂₃ -x ₁ y ₂ H ₃₁ -y ₁ y ₂ H ₃₂ -y ₂ H ₃₃ = 0 (23)

From this it is deduced that:

Since the last digit H ₃₃ of the perspective transformation model is usually set to 1 by us, the transformation model consists of eight parameters, and at least four pairs of non-collinear matching points are needed to obtain the values of the parameters of the perspective transformation model. However, no matter in the actual calculation process or in objective conditions, errors must exist. First of all, in actual calculation, due to noise interference, positioning error of feature point detection, selected model error, mismatch and other factors, if only four pairs of non-collinear matching point pairs are used to obtain eight parameters, the error will be large, and it is impossible To ensure the registration accuracy, the RANSAC algorithm can just solve the above problems very well.

As a selection mechanism, the RANSAC algorithm not only guarantees the number and quality of extracted matching point pairs as much as possible, but also accurately estimates the transformation model. The steps to obtain the parameters of the transformation matrix by the RANSAC algorithm are as follows:

(1) Given two color remote sensing images image1 and image2, find the pair of feature points m _i '←→m _i after two-way matching and purification, i=1,2,...n;

(2) From the set of matching points, randomly select 4 pairs of matching points, and use the coordinate values of these 4 pairs of matching points to obtain the transformation matrix H;

(3) Based on the Euclidean distance, find a point pair in the matching point pair set that meets the condition of d(Hm _i ,m _i ')<t, where Hm _i represents the matrix transformation mapping of the feature points in image2 to image1 The position coordinates in , m _i ' represents the position coordinates of _mi corresponding feature points in image1, d(Hm _i , m _i ') represents the Euclidean distance between two coordinate vectors, and the qualified point pairs are used as the final inner point, and record the number of interior points that satisfy the H constraint, and select 4 for t;

(4) Repeat (2) and (3) two steps n times (100≤n≤200), record the number of interior points each time;

(5) Select H _best corresponding to the largest number of interior points, find all matching point pairs that meet the condition of d(H _best m _i ,m _i ')<t, and use them as the final interior points, that is, the correct matching point pairs, Wrong matching point pairs that do not meet the condition of d(H _best m _i ,m _i ')<t are outliers and should be eliminated;

(6) According to the random sampling consensus algorithm to obtain N matching point pairs, we denote the matrix formed by 2N equations including matching point pairs and the corresponding transformation matrix H as A, and perform SVD singular value decomposition on it to obtain A= UDV ^T , U and V are the singular vectors of A, D is the singular values on the diagonal in descending order, and the last column of V is reconstructed into a 3*3 matrix, which is the parameter estimation of the perspective transformation matrix. This The parameters of the matrix satisfy the minimum error in the sense of least squares;

Step 5. Bicubic interpolation to achieve registration

After obtaining the transformation matrix, it is necessary to perform mapping to achieve registration. During the mapping process, since the pixels are stored in integer coordinate values, which are discontinuous values, this will lead to an inevitable situation: that is The original point (x ₁ , y ₁ ) (coordinates are all integers) on the integer grid does not fall on the grid point (x ₂ , y ₂ ) after mapping, so the role of interpolation is particularly important.

The bicubic interpolation method uses the gray information of sixteen pixels within the 4*4 neighborhood of the corresponding point in the image to be registered to estimate the gray value of the point to be interpolated, not only considering the four directly adjacent The influence of the gray value of each point is also considered, and the influence of the change rate of the gray value of each adjacent point is also considered. This method can be used if there are sixteen pixel points within a 4*4 neighborhood near the corresponding point in the image to be registered for the point to be interpolated.

If the coordinates of any point on the reference image are (x, y), its gray value is g(x, y), and the coordinates of the corresponding point on the image to be registered are (u, v), and its gray value is f(u,v), [ ] means taking an integer, and the matrix composed of 16 pixels in the 4*4 neighborhood of the corresponding point in the image to be registered is B:

The function f(u+i,v+j) (i=-1,0,1,2; j=-1,0,1,2) in the elements in the matrix B is defined as the corresponding point on the image to be registered The gray value when the coordinates are (u+i, v+j).

Then the calculation formula of the gray value of the point to be interpolated in the reference image is as follows:

g(x,y)=f(u,v)=ABC (26)

in:

s(w) is the basis function of bicubic interpolation, which is an approximation to sin(x*π)/x. Through bicubic interpolation, the image2 image is interpolated into the reference image image1 to achieve the final high-precision registration.

The bicubic interpolation method uses the grayscale information of the pixels in the larger neighborhood, the interpolation accuracy is further improved, and the edge details of the image are also well maintained, and the visual effect is more ideal. Basically, the interpolation accuracy is the highest among the existing interpolation algorithms algorithm.

The following example uses the original SURF algorithm and a robust SURF algorithm to verify the accuracy of the same two UAV color remote sensing images image1 and image2. and Figure 3). Between image1 and image2 is actually an affine transformation that has been scaled by 2 times, rotated by 10 degrees, offset by 239 pixels horizontally, and offset by 28 pixels vertically.

Example:

Firstly, according to the geometric transformation relationship between the color remote sensing images image1 and image2 mentioned above, we can easily obtain the theoretical value H _ture . Secondly, we use the original SURF algorithm to register the color remote sensing images image1 and image2 to obtain the original H _surf . Finally, a robust SURF algorithm is used to register the color remote sensing images image1 and image2 to obtain H _rob .

Without loss of generality, eight groups of UAV color remote sensing images were selected for experiments respectively. Here, eight groups of experiments are carried out to obtain the average value, and the statistics of the parameters of the transformation matrix are shown in Table 1. Table 1 is the parameters of the affine transformation matrix Comparing the table, it is found that the average value of the maximum deviation of the SURF algorithm is 0.218, which is higher than the average value of the maximum deviation of the robust algorithm, which is 0.134. The degree of proximity of the six parameters, the robust algorithm is significantly better than the SURF algorithm.

Table 1 Affine transformation matrix parameter comparison table

Claims (4)

1. A method for registration of robust SURF unmanned aerial vehicle color remote sensing images, is characterized in that the steps are as follows: Step 1: Perform grayscale processing on the acquired UAV color remote sensing reference image image1 and image2 to be registered to become a grayscale image, and perform Gaussian smoothing on the grayscale image to complete the scale space construction, and use the feature detection operator Hessian matrix The determinant of the feature points is extracted, and all the feature points of image1 and image2 are recorded as the target set {p _i }, i=1,2,...,n and the reference set {q _j },j=1,2, ...,m; Step 2: For the obtained target set {p _i } and reference set {q _j }, a square area with a side length of 20s is determined centering on each feature point, where s represents the sampling interval in different scale spaces, and the surrounding Neighborhood information is divided into 4×4 sub-regions, and each small region is divided into 5×5 sampling points. Finally, Harr wavelet is used to calculate the response and modulus dx,|dx|,dy, |dy|, and counting the total response Σdx, ∑|dx|, ∑dy, ∑|dy| of 5×5 sampling points, the 64-dimensional SURF feature descriptor of the feature point can be obtained: V=(∑dx,∑|dx|,∑dy,∑|dy|) Then the SURF descriptor of the feature point is defined in the horizontal or vertical direction as: V _gray ＝(V ₁ ,V ₂ ,V ₃ ,V ₄ ,...,V ₁₆ ) Among them, V _i , _{i =} 1, 2, 3,..., 15, 16 represent the description vectors corresponding to the 16 sub-regions; from this we can get all the 64-dimensional SURF feature descriptor for feature points; Step 3: For the feature point with coordinates (x, y), it is expressed as R _{(x, y)} , G _{(x, y)} , B _{(x, y)} in the RGB color space, and it is converted into the color space, and the conversion To the YIQ color space is Y _(x,y) , I _(x,y) , Q _(x,y) ; The transformation from RGB color space to YIQ color space is expressed as: $\left[\begin{array}{c}\mathrm{<span\; class="notranslate">\; Y</span>}\\ \mathrm{<span\; class="notranslate">\; I</span>}\\ \mathrm{<span\; class="notranslate">\; Q</span>}\end{array}\right]<span\; class="notranslate">\; =</span>\left[\begin{array}{ccc}\mathrm{<span\; class="notranslate">\; 0.229</span>}& \mathrm{<span\; class="notranslate">\; 0.587</span>}& \mathrm{<span\; class="notranslate">\; 0.114</span>}\\ \mathrm{<span\; class="notranslate">\; 0.596</span>}& <span\; class="notranslate">\; -</span>\mathrm{<span\; class="notranslate">\; 0.274</span>}& <span\; class="notranslate">\; -</span>\mathrm{<span\; class="notranslate">\; 0.322</span>}\\ \mathrm{<span\; class="notranslate">\; 0.211</span>}& <span\; class="notranslate">\; -</span>\mathrm{<span\; class="notranslate">\; 0.523</span>}& \mathrm{<span\; class="notranslate">\; 0.312</span>}\end{array}\right]\left[\begin{array}{c}\mathrm{<span\; class="notranslate">\; R</span>}\\ \mathrm{<span\; class="notranslate">\; G</span>}\\ \mathrm{<span\; class="notranslate">\; B</span>}\end{array}\right]$ Adding it to the original algorithm descriptor, the new color descriptor is as follows: ${\mathrm{<span\; class="notranslate">\; V</span>}}_{\mathrm{<span\; class="notranslate">\; c</span>}\mathrm{<span\; class="notranslate">\; o</span>}\mathrm{<span\; class="notranslate">\; l</span>}\mathrm{<span\; class="notranslate">\; o</span>}\mathrm{<span\; class="notranslate">\; r</span>}}<span\; class="notranslate">\; =</span><span\; class="notranslate">\; (</span>{\mathrm{<span\; class="notranslate">\; i</span>}}_{\mathrm{<span\; class="notranslate">\; 1</span>}}<span\; class="notranslate">\; ,</span>{\mathrm{<span\; class="notranslate">\; i</span>}}_{\mathrm{<span\; class="notranslate">\; 2</span>}}<span\; class="notranslate">\; ,</span>{\mathrm{<span\; class="notranslate">\; i</span>}}_{\mathrm{<span\; class="notranslate">\; 3</span>}}<span\; class="notranslate">\; ,</span>{\mathrm{<span\; class="notranslate">\; i</span>}}_{\mathrm{<span\; class="notranslate">\; 4</span>}}<span\; class="notranslate">\; ,</span><span\; class="notranslate">\; ...</span><span\; class="notranslate">\; ,</span>{\mathrm{<span\; class="notranslate">\; i</span>}}_{\mathrm{<span\; class="notranslate">\; 64</span>}}<span\; class="notranslate">\; ,</span>{\mathrm{<span\; class="notranslate">\; Y</span>}}_{<span\; class="notranslate">\; (</span>{\mathrm{<span\; class="notranslate">\; x</span>}}_{\mathrm{<span\; class="notranslate">\; 0</span>}}<span\; class="notranslate">\; ,</span>{\mathrm{<span\; class="notranslate">\; the\; y</span>}}_{\mathrm{<span\; class="notranslate">\; 0</span>}}<span\; class="notranslate">\; )</span>}<span\; class="notranslate">\; ,</span>{\mathrm{<span\; class="notranslate">\; I</span>}}_{<span\; class="notranslate">\; (</span>{\mathrm{<span\; class="notranslate">\; x</span>}}_{\mathrm{<span\; class="notranslate">\; 0</span>}}<span\; class="notranslate">\; ,</span>{\mathrm{<span\; class="notranslate">\; the\; y</span>}}_{\mathrm{<span\; class="notranslate">\; 0</span>}}<span\; class="notranslate">\; )</span>}<span\; class="notranslate">\; ,</span>{\mathrm{<span\; class="notranslate">\; Q</span>}}_{<span\; class="notranslate">\; (</span>{\mathrm{<span\; class="notranslate">\; x</span>}}_{\mathrm{<span\; class="notranslate">\; 0</span>}}<span\; class="notranslate">\; ,</span>{\mathrm{<span\; class="notranslate">\; the\; y</span>}}_{\mathrm{<span\; class="notranslate">\; 0</span>}}<span\; class="notranslate">\; )</span>}<span\; class="notranslate">\; )</span>$ Normalize the above formula: ${\stackrel{<span\; class="notranslate">\; \&OverBar;</span>}{\mathrm{<span\; class="notranslate">\; V</span>}}}_{\mathrm{<span\; class="notranslate">\; c</span>}\mathrm{<span\; class="notranslate">\; o</span>}\mathrm{<span\; class="notranslate">\; l</span>}\mathrm{<span\; class="notranslate">\; o</span>}\mathrm{<span\; class="notranslate">\; r</span>}}<span\; class="notranslate">\; =</span><span\; class="notranslate">\; (</span>\frac{{\mathrm{<span\; class="notranslate">\; i</span>}}_{\mathrm{<span\; class="notranslate">\; 1</span>}}}{<span\; class="notranslate">\; |</span>{\mathrm{<span\; class="notranslate">\; V</span>}}_{\mathrm{<span\; class="notranslate">\; c</span>}\mathrm{<span\; class="notranslate">\; o</span>}\mathrm{<span\; class="notranslate">\; l</span>}\mathrm{<span\; class="notranslate">\; o</span>}\mathrm{<span\; class="notranslate">\; r</span>}}<span\; class="notranslate">\; |</span>}<span\; class="notranslate">\; ,</span>\frac{{\mathrm{<span\; class="notranslate">\; i</span>}}_{\mathrm{<span\; class="notranslate">\; 2</span>}}}{<span\; class="notranslate">\; |</span>{\mathrm{<span\; class="notranslate">\; V</span>}}_{\mathrm{<span\; class="notranslate">\; c</span>}\mathrm{<span\; class="notranslate">\; o</span>}\mathrm{<span\; class="notranslate">\; l</span>}\mathrm{<span\; class="notranslate">\; o</span>}\mathrm{<span\; class="notranslate">\; r</span>}}<span\; class="notranslate">\; |</span>}<span\; class="notranslate">\; ,</span>\frac{{\mathrm{<span\; class="notranslate">\; i</span>}}_{\mathrm{<span\; class="notranslate">\; 3</span>}}}{<span\; class="notranslate">\; |</span>{\mathrm{<span\; class="notranslate">\; V</span>}}_{\mathrm{<span\; class="notranslate">\; c</span>}\mathrm{<span\; class="notranslate">\; o</span>}\mathrm{<span\; class="notranslate">\; l</span>}\mathrm{<span\; class="notranslate">\; o</span>}\mathrm{<span\; class="notranslate">\; r</span>}}<span\; class="notranslate">\; |</span>}<span\; class="notranslate">\; ,</span>\frac{{\mathrm{<span\; class="notranslate">\; i</span>}}_{\mathrm{<span\; class="notranslate">\; 4</span>}}}{<span\; class="notranslate">\; |</span>{\mathrm{<span\; class="notranslate">\; V</span>}}_{\mathrm{<span\; class="notranslate">\; c</span>}\mathrm{<span\; class="notranslate">\; o</span>}\mathrm{<span\; class="notranslate">\; l</span>}\mathrm{<span\; class="notranslate">\; o</span>}\mathrm{<span\; class="notranslate">\; r</span>}}<span\; class="notranslate">\; |</span>}<span\; class="notranslate">\; ,</span><span\; class="notranslate">\; ...</span><span\; class="notranslate">\; ,</span>\frac{{\mathrm{<span\; class="notranslate">\; i</span>}}_{\mathrm{<span\; class="notranslate">\; 64</span>}}}{<span\; class="notranslate">\; |</span>{\mathrm{<span\; class="notranslate">\; V</span>}}_{\mathrm{<span\; class="notranslate">\; c</span>}\mathrm{<span\; class="notranslate">\; o</span>}\mathrm{<span\; class="notranslate">\; l</span>}\mathrm{<span\; class="notranslate">\; o</span>}\mathrm{<span\; class="notranslate">\; r</span>}}<span\; class="notranslate">\; |</span>}<span\; class="notranslate">\; ,</span>\frac{{\mathrm{<span\; class="notranslate">\; Y</span>}}_{<span\; class="notranslate">\; (</span>{\mathrm{<span\; class="notranslate">\; x</span>}}_{\mathrm{<span\; class="notranslate">\; 0</span>}}<span\; class="notranslate">\; ,</span>{\mathrm{<span\; class="notranslate">\; the\; y</span>}}_{\mathrm{<span\; class="notranslate">\; 0</span>}}<span\; class="notranslate">\; )</span>}}{<span\; class="notranslate">\; |</span>{\mathrm{<span\; class="notranslate">\; V</span>}}_{\mathrm{<span\; class="notranslate">\; c</span>}\mathrm{<span\; class="notranslate">\; o</span>}\mathrm{<span\; class="notranslate">\; l</span>}\mathrm{<span\; class="notranslate">\; o</span>}\mathrm{<span\; class="notranslate">\; r</span>}}<span\; class="notranslate">\; |</span>}<span\; class="notranslate">\; ,</span>\frac{{\mathrm{<span\; class="notranslate">\; I</span>}}_{<span\; class="notranslate">\; (</span>{\mathrm{<span\; class="notranslate">\; x</span>}}_{\mathrm{<span\; class="notranslate">\; 0</span>}}<span\; class="notranslate">\; ,</span>{\mathrm{<span\; class="notranslate">\; the\; y</span>}}_{\mathrm{<span\; class="notranslate">\; 0</span>}}<span\; class="notranslate">\; )</span>}}{<span\; class="notranslate">\; |</span>{\mathrm{<span\; class="notranslate">\; V</span>}}_{\mathrm{<span\; class="notranslate">\; c</span>}\mathrm{<span\; class="notranslate">\; o</span>}\mathrm{<span\; class="notranslate">\; l</span>}\mathrm{<span\; class="notranslate">\; o</span>}\mathrm{<span\; class="notranslate">\; r</span>}}<span\; class="notranslate">\; |</span>}<span\; class="notranslate">\; ,</span>\frac{{\mathrm{<span\; class="notranslate">\; Q</span>}}_{<span\; class="notranslate">\; (</span>{\mathrm{<span\; class="notranslate">\; x</span>}}_{\mathrm{<span\; class="notranslate">\; 0</span>}}<span\; class="notranslate">\; ,</span>{\mathrm{<span\; class="notranslate">\; the\; y</span>}}_{\mathrm{<span\; class="notranslate">\; 0</span>}}<span\; class="notranslate">\; )</span>}}{<span\; class="notranslate">\; |</span>{\mathrm{<span\; class="notranslate">\; V</span>}}_{\mathrm{<span\; class="notranslate">\; c</span>}\mathrm{<span\; class="notranslate">\; o</span>}\mathrm{<span\; class="notranslate">\; l</span>}\mathrm{<span\; class="notranslate">\; o</span>}\mathrm{<span\; class="notranslate">\; r</span>}}<span\; class="notranslate">\; |</span>}<span\; class="notranslate">\; )</span>$ $<span\; class="notranslate">\; |</span>{\mathrm{<span\; class="notranslate">\; V</span>}}_{\mathrm{<span\; class="notranslate">\; c</span>}\mathrm{<span\; class="notranslate">\; o</span>}\mathrm{<span\; class="notranslate">\; l</span>}\mathrm{<span\; class="notranslate">\; o</span>}\mathrm{<span\; class="notranslate">\; r</span>}}<span\; class="notranslate">\; |</span><span\; class="notranslate">\; =</span>\sqrt{\underset{\mathrm{<span\; class="notranslate">\; k</span>}<span\; class="notranslate">\; =</span>\mathrm{<span\; class="notranslate">\; 1</span>}}{\overset{\mathrm{<span\; class="notranslate">\; 64</span>}}{<span\; class="notranslate">\; \&Sigma;</span>}}{{\mathrm{<span\; class="notranslate">\; i</span>}}_{\mathrm{<span\; class="notranslate">\; k</span>}}}^{\mathrm{<span\; class="notranslate">\; 2</span>}}<span\; class="notranslate">\; +</span>{{\mathrm{<span\; class="notranslate">\; Y</span>}}^{\mathrm{<span\; class="notranslate">\; 2</span>}}}_{<span\; class="notranslate">\; (</span>{\mathrm{<span\; class="notranslate">\; x</span>}}_{\mathrm{<span\; class="notranslate">\; 0</span>}}<span\; class="notranslate">\; ,</span>{\mathrm{<span\; class="notranslate">\; the\; y</span>}}_{\mathrm{<span\; class="notranslate">\; 0</span>}}<span\; class="notranslate">\; )</span>}<span\; class="notranslate">\; +</span>{{\mathrm{<span\; class="notranslate">\; I</span>}}^{\mathrm{<span\; class="notranslate">\; 2</span>}}}_{<span\; class="notranslate">\; (</span>{\mathrm{<span\; class="notranslate">\; x</span>}}_{\mathrm{<span\; class="notranslate">\; 0</span>}}<span\; class="notranslate">\; ,</span>{\mathrm{<span\; class="notranslate">\; the\; y</span>}}_{\mathrm{<span\; class="notranslate">\; 0</span>}}<span\; class="notranslate">\; )</span>}<span\; class="notranslate">\; +</span>{{\mathrm{<span\; class="notranslate">\; Q</span>}}^{\mathrm{<span\; class="notranslate">\; 2</span>}}}_{<span\; class="notranslate">\; (</span>{\mathrm{<span\; class="notranslate">\; x</span>}}_{\mathrm{<span\; class="notranslate">\; 0</span>}}<span\; class="notranslate">\; ,</span>{\mathrm{<span\; class="notranslate">\; the\; y</span>}}_{\mathrm{<span\; class="notranslate">\; 0</span>}}<span\; class="notranslate">\; )</span>}}$ where _{i k} , _k =1,2,3,...,63,64 represent the 64-dimensional description vector of feature points in the classic SURF algorithm; 67-dimensional description vector of each feature point; Step 4: Two-way matching of feature points Step 4a: use the Euclidean distance formula to calculate the nearest Euclidean distance d _nt and the next closest Euclidean distance d _snt in the forward direction of the feature point p1 in the reference image image1; Use the Euclidean distance formula to calculate the feature points of the feature point p1 in the reference image image1 in the forward direction of the nearest Euclidean distance d1 _nt and the second closest Euclidean distance d1 _snt in the image2, respectively, and the feature points are p2 _nt and p2 _snt ; Step 4b: Calculate the ratio T ₁ of the closest distance d1 _nt and the second closest distance d1 _snt , namely ${\mathrm{<span\; class="notranslate">\; T</span>}}_{\mathrm{<span\; class="notranslate">\; 1</span>}}<span\; class="notranslate">\; =</span>\frac{\mathrm{<span\; class="notranslate">\; d</span>}{\mathrm{<span\; class="notranslate">\; 1</span>}}_{\mathrm{<span\; class="notranslate">\; no</span>}\mathrm{<span\; class="notranslate">\; t</span>}}}{\mathrm{<span\; class="notranslate">\; d</span>}{\mathrm{<span\; class="notranslate">\; 1</span>}}_{\mathrm{<span\; class="notranslate">\; the\; s</span>}\mathrm{<span\; class="notranslate">\; no</span>}\mathrm{<span\; class="notranslate">\; t</span>}}}$ Compare T ₁ with the judgment threshold T: if T ₁ <T, proceed to step 4c, otherwise, jump out automatically; Step 4c: For the feature point corresponding to p2 _nt , use the Euclidean distance formula to calculate the feature points of the reverse closest distance d2 _nt and the next closest distance d2 _snt in image1 are p1 _nt and p1 _snt respectively; Step 4d: Calculate the ratio T ₂ of the closest distance to the next closest distance ${\mathrm{<span\; class="notranslate">\; T</span>}}_{\mathrm{<span\; class="notranslate">\; 2</span>}}<span\; class="notranslate">\; =</span>\frac{\mathrm{<span\; class="notranslate">\; d</span>}{\mathrm{<span\; class="notranslate">\; 2</span>}}_{\mathrm{<span\; class="notranslate">\; no</span>}\mathrm{<span\; class="notranslate">\; t</span>}}}{\mathrm{<span\; class="notranslate">\; d</span>}{\mathrm{<span\; class="notranslate">\; 2</span>}}_{\mathrm{<span\; class="notranslate">\; the\; s</span>}\mathrm{<span\; class="notranslate">\; no</span>}\mathrm{<span\; class="notranslate">\; t</span>}}}$ Compare T ₁ with the judgment threshold T: if T ₂ <T, and the feature point p1 _nt is looking for is the same feature point as the original feature point p1, it is considered to be a correct matching point pair, otherwise return to step 4a for the reference image Other feature points in image1 are used to judge new feature points; Traverse all the feature points in the reference image target set in turn, and finally achieve two-way matching and purification; Step 5: Use the RANSAC algorithm to obtain the parameters of the transformation matrix Step 5a: Find the feature point pair m _i '←→m _i after the two-way matching and purification after step 4 is completed, i=1,2,...n; Step 5b: Randomly select 4 pairs of matching points from the set of matching points, and use the coordinate values of these 4 pairs of matching points to obtain the transformation matrix H; Step 5c: Based on the Euclidean distance, find the point pair in the matching point pair set that meets the condition of d(H _i m _i ,m _i ')<t, where H _i m _i represents the matrix of the feature points in image2 The transformation is mapped to the position coordinates in image1, m _i 'indicates the position coordinates of _mi in image1 corresponding to the feature point, d(H _i m _i ,m _i ') indicates the Euclidean distance between two coordinate vectors, and will meet the conditions The point pair is used as the final interior point, and the number of interior points satisfying the H _i constraint is recorded; Step 5d: Repeat steps 5b and 5c n times, and record the number of interior points each time; Step 5e: Select H _best corresponding to the largest number of interior points, find all matching point pairs that meet the condition of d(H _best m _i ,m _i ')<t, and use them as the final interior points, that is, the correct matching point pairs, Wrong matching point pairs that do not meet the condition of d(H _best m _i ,m _i ')<t are outliers and should be eliminated; Step 5f: According to the random sampling consensus algorithm, N matching point pairs are obtained, and the 2N matrices composed of x ₁ , y ₁ , x ₂ , y ₂ are marked, denoted as A, and the SVD singular value is performed on it Decompose A=UDV ^T , U and V are the singular vectors of A, D is the singular values on the diagonal in descending order, and the last column of V is reconstructed into a 3*3 matrix, which is the required perspective transformation matrix Parameter Estimation; Step 6: Bicubic interpolation for registration If the coordinates of any point on the reference image image1 are (x, y), its gray value is g(x, y), and the corresponding point coordinates of this point on the image to be registered are (u, v), and its gray value The value is f(u,v), and [ ] means an integer, and the matrix composed of 16 pixels in the 4*4 neighborhood of the corresponding point in the image to be registered is B: $\mathrm{<span\; class="notranslate">\; B</span>}<span\; class="notranslate">\; =</span>\left[\begin{array}{cccc}\mathrm{<span\; class="notranslate">\; f</span>}<span\; class="notranslate">\; (</span><span\; class="notranslate">\; \&lsqb;</span>\mathrm{<span\; class="notranslate">\; u</span>}<span\; class="notranslate">\; \&rsqb;</span><span\; class="notranslate">\; -</span>\mathrm{<span\; class="notranslate">\; 1</span>}<span\; class="notranslate">\; ,</span><span\; class="notranslate">\; \&lsqb;</span>\mathrm{<span\; class="notranslate">\; v</span>}<span\; class="notranslate">\; \&rsqb;</span><span\; class="notranslate">\; -</span>\mathrm{<span\; class="notranslate">\; 1</span>}<span\; class="notranslate">\; )</span>& \mathrm{<span\; class="notranslate">\; f</span>}<span\; class="notranslate">\; (</span><span\; class="notranslate">\; \&lsqb;</span>\mathrm{<span\; class="notranslate">\; u</span>}<span\; class="notranslate">\; \&rsqb;</span><span\; class="notranslate">\; -</span>\mathrm{<span\; class="notranslate">\; 1</span>}<span\; class="notranslate">\; ,</span><span\; class="notranslate">\; \&lsqb;</span>\mathrm{<span\; class="notranslate">\; v</span>}<span\; class="notranslate">\; \&rsqb;</span><span\; class="notranslate">\; )</span>& \mathrm{<span\; class="notranslate">\; f</span>}<span\; class="notranslate">\; (</span><span\; class="notranslate">\; \&lsqb;</span>\mathrm{<span\; class="notranslate">\; u</span>}<span\; class="notranslate">\; \&rsqb;</span><span\; class="notranslate">\; -</span>\mathrm{<span\; class="notranslate">\; 1</span>}<span\; class="notranslate">\; ,</span><span\; class="notranslate">\; \&lsqb;</span>\mathrm{<span\; class="notranslate">\; v</span>}<span\; class="notranslate">\; \&rsqb;</span><span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; 1</span>}<span\; class="notranslate">\; )</span>& \mathrm{<span\; class="notranslate">\; f</span>}<span\; class="notranslate">\; (</span><span\; class="notranslate">\; \&lsqb;</span>\mathrm{<span\; class="notranslate">\; u</span>}<span\; class="notranslate">\; \&rsqb;</span><span\; class="notranslate">\; -</span>\mathrm{<span\; class="notranslate">\; 1</span>}<span\; class="notranslate">\; ,</span><span\; class="notranslate">\; \&lsqb;</span>\mathrm{<span\; class="notranslate">\; v</span>}<span\; class="notranslate">\; \&rsqb;</span><span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; 2</span>}<span\; class="notranslate">\; )</span>\\ \mathrm{<span\; class="notranslate">\; f</span>}<span\; class="notranslate">\; (</span><span\; class="notranslate">\; \&lsqb;</span>\mathrm{<span\; class="notranslate">\; u</span>}<span\; class="notranslate">\; \&rsqb;</span><span\; class="notranslate">\; ,</span><span\; class="notranslate">\; \&lsqb;</span>\mathrm{<span\; class="notranslate">\; v</span>}<span\; class="notranslate">\; \&rsqb;</span><span\; class="notranslate">\; -</span>\mathrm{<span\; class="notranslate">\; 1</span>}<span\; class="notranslate">\; )</span>& \mathrm{<span\; class="notranslate">\; f</span>}<span\; class="notranslate">\; (</span><span\; class="notranslate">\; \&lsqb;</span>\mathrm{<span\; class="notranslate">\; u</span>}<span\; class="notranslate">\; \&rsqb;</span><span\; class="notranslate">\; ,</span><span\; class="notranslate">\; \&lsqb;</span>\mathrm{<span\; class="notranslate">\; v</span>}<span\; class="notranslate">\; \&rsqb;</span><span\; class="notranslate">\; )</span>& \mathrm{<span\; class="notranslate">\; f</span>}<span\; class="notranslate">\; (</span><span\; class="notranslate">\; \&lsqb;</span>\mathrm{<span\; class="notranslate">\; u</span>}<span\; class="notranslate">\; \&rsqb;</span><span\; class="notranslate">\; ,</span><span\; class="notranslate">\; \&lsqb;</span>\mathrm{<span\; class="notranslate">\; v</span>}<span\; class="notranslate">\; \&rsqb;</span><span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; 1</span>}<span\; class="notranslate">\; )</span>& \mathrm{<span\; class="notranslate">\; f</span>}<span\; class="notranslate">\; (</span><span\; class="notranslate">\; \&lsqb;</span>\mathrm{<span\; class="notranslate">\; u</span>}<span\; class="notranslate">\; \&rsqb;</span><span\; class="notranslate">\; ,</span><span\; class="notranslate">\; \&lsqb;</span>\mathrm{<span\; class="notranslate">\; v</span>}<span\; class="notranslate">\; \&rsqb;</span><span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; 2</span>}<span\; class="notranslate">\; )</span>\\ \mathrm{<span\; class="notranslate">\; f</span>}<span\; class="notranslate">\; (</span><span\; class="notranslate">\; \&lsqb;</span>\mathrm{<span\; class="notranslate">\; u</span>}<span\; class="notranslate">\; \&rsqb;</span><span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; 1</span>}<span\; class="notranslate">\; ,</span><span\; class="notranslate">\; \&lsqb;</span>\mathrm{<span\; class="notranslate">\; v</span>}<span\; class="notranslate">\; \&rsqb;</span><span\; class="notranslate">\; -</span>\mathrm{<span\; class="notranslate">\; 1</span>}<span\; class="notranslate">\; )</span>& \mathrm{<span\; class="notranslate">\; f</span>}<span\; class="notranslate">\; (</span><span\; class="notranslate">\; \&lsqb;</span>\mathrm{<span\; class="notranslate">\; u</span>}<span\; class="notranslate">\; \&rsqb;</span><span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; 1</span>}<span\; class="notranslate">\; ,</span><span\; class="notranslate">\; \&lsqb;</span>\mathrm{<span\; class="notranslate">\; v</span>}<span\; class="notranslate">\; \&rsqb;</span><span\; class="notranslate">\; )</span>& \mathrm{<span\; class="notranslate">\; f</span>}<span\; class="notranslate">\; (</span><span\; class="notranslate">\; \&lsqb;</span>\mathrm{<span\; class="notranslate">\; u</span>}<span\; class="notranslate">\; \&rsqb;</span><span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; 1</span>}<span\; class="notranslate">\; ,</span><span\; class="notranslate">\; \&lsqb;</span>\mathrm{<span\; class="notranslate">\; v</span>}<span\; class="notranslate">\; \&rsqb;</span><span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; 1</span>}<span\; class="notranslate">\; )</span>& \mathrm{<span\; class="notranslate">\; f</span>}<span\; class="notranslate">\; (</span><span\; class="notranslate">\; \&lsqb;</span>\mathrm{<span\; class="notranslate">\; u</span>}<span\; class="notranslate">\; \&rsqb;</span><span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; 1</span>}<span\; class="notranslate">\; ,</span><span\; class="notranslate">\; \&lsqb;</span>\mathrm{<span\; class="notranslate">\; v</span>}<span\; class="notranslate">\; \&rsqb;</span><span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; 2</span>}<span\; class="notranslate">\; )</span>\\ \mathrm{<span\; class="notranslate">\; f</span>}<span\; class="notranslate">\; (</span><span\; class="notranslate">\; \&lsqb;</span>\mathrm{<span\; class="notranslate">\; u</span>}<span\; class="notranslate">\; \&rsqb;</span><span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; 2</span>}<span\; class="notranslate">\; ,</span><span\; class="notranslate">\; \&lsqb;</span>\mathrm{<span\; class="notranslate">\; v</span>}<span\; class="notranslate">\; \&rsqb;</span><span\; class="notranslate">\; -</span>\mathrm{<span\; class="notranslate">\; 1</span>}<span\; class="notranslate">\; )</span>& \mathrm{<span\; class="notranslate">\; f</span>}<span\; class="notranslate">\; (</span><span\; class="notranslate">\; \&lsqb;</span>\mathrm{<span\; class="notranslate">\; u</span>}<span\; class="notranslate">\; \&rsqb;</span><span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; 2</span>}<span\; class="notranslate">\; ,</span><span\; class="notranslate">\; \&lsqb;</span>\mathrm{<span\; class="notranslate">\; v</span>}<span\; class="notranslate">\; \&rsqb;</span><span\; class="notranslate">\; )</span>& \mathrm{<span\; class="notranslate">\; f</span>}<span\; class="notranslate">\; (</span><span\; class="notranslate">\; \&lsqb;</span>\mathrm{<span\; class="notranslate">\; u</span>}<span\; class="notranslate">\; \&rsqb;</span><span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; 2</span>}<span\; class="notranslate">\; ,</span><span\; class="notranslate">\; \&lsqb;</span>\mathrm{<span\; class="notranslate">\; v</span>}<span\; class="notranslate">\; \&rsqb;</span><span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; 1</span>}<span\; class="notranslate">\; )</span>& \mathrm{<span\; class="notranslate">\; f</span>}<span\; class="notranslate">\; (</span><span\; class="notranslate">\; \&lsqb;</span>\mathrm{<span\; class="notranslate">\; u</span>}<span\; class="notranslate">\; \&rsqb;</span><span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; 2</span>}<span\; class="notranslate">\; ,</span><span\; class="notranslate">\; \&lsqb;</span>\mathrm{<span\; class="notranslate">\; v</span>}<span\; class="notranslate">\; \&rsqb;</span><span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; 2</span>}<span\; class="notranslate">\; )</span>\end{array}\right]$ The function f(u+i,v+j) in the elements in the matrix B, i=-1,0,1,2; j=-1,0,1,2, defined as the corresponding point on the image to be registered The gray value when the coordinates are (u+i,v+j); Then the calculation formula of the gray value of the point to be interpolated in the reference image is as follows: g(x,y)=f(u,v)=ABC in: $\left\{\begin{array}{c}\mathrm{<span\; class="notranslate">\; a</span>}<span\; class="notranslate">\; =</span>\mathrm{<span\; class="notranslate">\; u</span>}<span\; class="notranslate">\; -</span><span\; class="notranslate">\; \&lsqb;</span>\mathrm{<span\; class="notranslate">\; u</span>}<span\; class="notranslate">\; \&rsqb;</span>\\ \mathrm{<span\; class="notranslate">\; b</span>}<span\; class="notranslate">\; =</span>\mathrm{<span\; class="notranslate">\; v</span>}<span\; class="notranslate">\; -</span><span\; class="notranslate">\; \&lsqb;</span>\mathrm{<span\; class="notranslate">\; v</span>}<span\; class="notranslate">\; \&rsqb;</span>\\ \mathrm{<span\; class="notranslate">\; C</span>}<span\; class="notranslate">\; =</span>{\left[\begin{array}{cccc}\mathrm{<span\; class="notranslate">\; the\; s</span>}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; 1</span>}<span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; a</span>}<span\; class="notranslate">\; )</span>& \mathrm{<span\; class="notranslate">\; the\; s</span>}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; a</span>}<span\; class="notranslate">\; )</span>& \mathrm{<span\; class="notranslate">\; the\; s</span>}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; 1</span>}<span\; class="notranslate">\; -</span>\mathrm{<span\; class="notranslate">\; a</span>}<span\; class="notranslate">\; )</span>& \mathrm{<span\; class="notranslate">\; the\; s</span>}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; 2</span>}<span\; class="notranslate">\; -</span>\mathrm{<span\; class="notranslate">\; a</span>}<span\; class="notranslate">\; )</span>\end{array}\right]}^{\mathrm{<span\; class="notranslate">\; T</span>}}\\ \mathrm{<span\; class="notranslate">\; A</span>}<span\; class="notranslate">\; =</span>\left[\begin{array}{cccc}\mathrm{<span\; class="notranslate">\; the\; s</span>}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; 1</span>}<span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; b</span>}<span\; class="notranslate">\; )</span>& \mathrm{<span\; class="notranslate">\; the\; s</span>}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; b</span>}<span\; class="notranslate">\; )</span>& \mathrm{<span\; class="notranslate">\; the\; s</span>}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; 1</span>}<span\; class="notranslate">\; -</span>\mathrm{<span\; class="notranslate">\; b</span>}<span\; class="notranslate">\; )</span>& \mathrm{<span\; class="notranslate">\; the\; s</span>}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; 2</span>}<span\; class="notranslate">\; -</span>\mathrm{<span\; class="notranslate">\; b</span>}<span\; class="notranslate">\; )</span>\end{array}\right]\end{array}\right.$ $\mathrm{<span\; class="notranslate">\; the\; s</span>}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; w</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; =</span>\left\{\begin{array}{cc}\mathrm{<span\; class="notranslate">\; 1</span>}<span\; class="notranslate">\; -</span>\mathrm{<span\; class="notranslate">\; 2</span>}<span\; class="notranslate">\; |</span>\mathrm{<span\; class="notranslate">\; w</span>}{<span\; class="notranslate">\; |</span>}^{\mathrm{<span\; class="notranslate">\; 2</span>}}<span\; class="notranslate">\; +</span><span\; class="notranslate">\; |</span>\mathrm{<span\; class="notranslate">\; w</span>}{<span\; class="notranslate">\; |</span>}^{\mathrm{<span\; class="notranslate">\; 3</span>}}<span\; class="notranslate">\; ,</span>& <span\; class="notranslate">\; |</span>\mathrm{<span\; class="notranslate">\; w</span>}<span\; class="notranslate">\; |</span><span\; class="notranslate">\; <</span>\mathrm{<span\; class="notranslate">\; 1</span>}\\ \mathrm{<span\; class="notranslate">\; 4</span>}<span\; class="notranslate">\; -</span>\mathrm{<span\; class="notranslate">\; 8</span>}<span\; class="notranslate">\; |</span>\mathrm{<span\; class="notranslate">\; w</span>}<span\; class="notranslate">\; |</span><span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; 5</span>}<span\; class="notranslate">\; |</span>\mathrm{<span\; class="notranslate">\; w</span>}{<span\; class="notranslate">\; |</span>}^{\mathrm{<span\; class="notranslate">\; 2</span>}}<span\; class="notranslate">\; -</span><span\; class="notranslate">\; |</span>\mathrm{<span\; class="notranslate">\; w</span>}{<span\; class="notranslate">\; |</span>}^{\mathrm{<span\; class="notranslate">\; 3</span>}}<span\; class="notranslate">\; ,</span>& \mathrm{<span\; class="notranslate">\; 1</span>}<span\; class="notranslate">\; \&le;</span><span\; class="notranslate">\; |</span>\mathrm{<span\; class="notranslate">\; w</span>}<span\; class="notranslate">\; |</span><span\; class="notranslate">\; <</span>\mathrm{<span\; class="notranslate">\; 2</span>}\\ \mathrm{<span\; class="notranslate">\; 0</span>}<span\; class="notranslate">\; ,</span>& <span\; class="notranslate">\; |</span>\mathrm{<span\; class="notranslate">\; w</span>}<span\; class="notranslate">\; |</span><span\; class="notranslate">\; \&Greater\; Equal;</span>\mathrm{<span\; class="notranslate">\; 2</span>}\end{array}\right.$ s(w) is the basis function of bicubic interpolation, which is the approximation to sin(x*π)/x. Through bicubic interpolation, the image2 image is interpolated into the reference image image1 to achieve the final high-precision registration. 2. A kind of robust SURF unmanned aerial vehicle color remote sensing image registration method according to claim 1, it is characterized in that the judgment threshold T in the described step 4b is 0.4-0.8. 3. a kind of robust SURF UAV color remote sensing image registration method according to claim 1, is characterized in that t in the described step 5c is selected 4. 4. A method for registration of robust SURF UAV color remote sensing images according to claim 1, characterized in that 100≤n≤200 in the step 5d.