CN107808362A - A kind of image split-joint method combined based on unmanned plane POS information with image SURF features - Google Patents

Abstract

The invention discloses an image splicing method based on the combination of UAV POS information and image SURF features, and relates to the fields of digital image processing, GIS, surveying and mapping, and other related fields. First, the image is geometrically corrected, and then the geographic coordinates of the four corners of the image are calculated. Based on the geographic coordinates of the first image, the positional relationship of the matching target points with the same name is obtained by extracting the SURF features of the overlapping areas of adjacent images, and then corrected in turn. The geographic coordinates of the following images, and finally adopt the fusion algorithm of adaptive gradual in and out, and obtain a panoramic image with good visual effect, and complete the good splicing of the images. This method combines the image feature extraction algorithm with the geographic coordinates of the image, which has a greater improvement in stitching efficiency and visual effects than the traditional feature extraction stitching algorithm, and the stitched image has geographic information, which has certain practical value .

Description

An image mosaic based on the combination of drone POS information and image SURF features method

technical field

The invention relates to a splicing method based on the combination of unmanned aerial vehicle POS (Position and Orientation System) information and image SURF features, and belongs to the fields of digital image processing, GIS, surveying and mapping, and other related fields.

Background technique

UAV (Unmanned Aerial Vehicle) has the characteristics of simple operation, quick response, flexible flight, and low cost. It has been widely used in disaster relief, military reconnaissance, maritime detection, environmental protection, and land dynamic monitoring. However, due to the limited viewing angle of UAV aerial photography sequence pictures, it is limited by the flying height of UAV and the parameters of the camera. In order to fully grasp and analyze the captured area and get more information about the target area, it is very urgent and necessary to quickly stitch images.

At present, there are two main methods for stitching aerial images of drones: one is a stitching method based on image features, and the other is a stitching method based on UAV POS information.

Feature-based image stitching mainly includes two main steps, namely image registration and image fusion. Among them, image registration is the core step. In the image registration algorithm based on feature detection, in 1988, CHris Harris proposed the Harris corner detection operator. The SIFT feature matching algorithm (scale invariant feature transform) is a classic algorithm. By looking for extreme points in the spatial scale, and extracting its position, scale, and rotation invariants, and the proposed features are relatively stable, but the disadvantage is that the number of extracted feature points is large, the amount of calculation is large, and it takes a long time. Bay proposed the SURF algorithm in 2006, which means to speed up the robustness feature. It is an improvement on the SIFT algorithm. The SURF algorithm uses the Hessian matrix to determine the candidate points and then performs non-maximum suppression, which reduces the computational complexity. In other words, the standard SURF operator is several times faster than the SIFT operator, and it has better robustness under multiple pictures, but the error of the feature-based image stitching method is easy to accumulate, and it is not easy to target features such as sea areas. In the case of extraction, traditional feature stitching methods are not applicable.

The POS-based information splicing method mainly uses the geographical coordinate relationship of aerial images for registration. The flight speed, flight altitude, and aircraft latitude and longitude information of the current exposure point can all be obtained through the POS system carried by the UAV. . This information records the attitude of the aircraft at each exposure time when it is flying. However, since the attitude of the drone is not fixed when it is taking aerial photography, two sets of POS data of adjacent exposure points are needed for aerotriangulation calculation. Thus, the latitude and longitude of the pixel points of the image are calculated, and the error varies with the establishment of the model.

This method combines the feature extraction algorithm with the geographic coordinates of the image, which has a greater improvement in stitching efficiency and visual effects than the traditional feature extraction stitching algorithm, and the stitched image has geographic information, which has certain practical value.

Contents of the invention

In view of the shortcomings of the above two methods, in order to improve the stitching accuracy and increase the adaptability to various aerial images, a method based on the combination of UAV POS information and traditional feature stitching is proposed. First, the image is geometrically corrected. Then calculate the geographic coordinates of the four corners of the image. Based on the geographic coordinates of the first image, extract the SURF features of the overlapping areas of adjacent images to obtain the positional relationship of the matching target points with the same name, thereby sequentially correcting the geographic coordinates of the subsequent images. Finally By adopting the fusion algorithm of adaptive gradual in and gradual out, a panoramic image with good visual effect is obtained, and the good splicing of the images is completed.

The technical scheme adopted in the present invention is a method of image splicing based on the combination of UAV POS information and image SURF features, the method comprising the following steps:

S1: Aerial image preprocessing:

During the aerial photography mission of the UAV, due to the influence of the aircraft's attitude, height, speed, and the rotation of the earth, the image is squeezed, twisted, stretched, and offset relative to the actual position of the ground target. The images are geometrically rectified to obtain remote sensing images based on the same reference projection plane.

S2: Calculation of overlapping area of aerial images:

Calculate the overlapping area of adjacent images according to the POS information, obtain the exposure time of the sequence image and the actual distance of the flight according to the POS information, decompose the speed in the positive direction according to the current flight angle of the aircraft, and calculate the overlapping area of adjacent captured images In order to reduce the calculation amount of splicing, and as the search area of the feature points is smaller, the probability of false matching is smaller, and the detection efficiency is improved.

S3: Geographical coordinate calculation of the image:

After the image is geometrically corrected, the latitude and longitude coordinates of the center point of the image are calculated according to the attitude angle of the current aircraft flight, and the ground resolution of the image is calculated according to the internal and external orientation elements of the camera, and the latitude and longitude coordinates of the four points of the image are calculated, and then the geographic coordinates are converted To the space Cartesian coordinate system, so as to project according to the space Cartesian coordinate system, and obtain the relative positional relationship between the images.

S4: Extract the SURF feature of the image:

The extraction of SURF features includes: constructing a Hessian matrix, generating all interest points, constructing a scale space, locating feature points, assigning the main direction of feature points, generating feature point descriptors, and matching feature points. Through these steps, the adjacent The SURF features of the image match point pairs.

S5: Use the extracted image SURF feature point pairs to correct the geographic coordinates:

Since the geographic coordinates of each image are already known, the geographic coordinates of the obtained point pairs with the same name should be the same point. If the calculated ones are different, use the first image as the reference, and correct the subsequent geographic coordinates in order to determine the position. registration, improving the accuracy of placement directly by geographic coordinates.

S6: Image fusion strategy:

For the image after correcting the coordinates, since there are certain traces between each image, some strategies need to be adopted to solve the problem of large color transition differences between stitching gaps, and the gradual in and out adaptive weight fusion method is used to image The transition between the gaps is smooth, so that the pictures look natural and the visual effect is good.

Compared with the prior art, the present invention has the following advantages: the method proposed in this paper combines the image feature extraction algorithm with the geographic coordinates of the image, and has a greater improvement in stitching efficiency and visual effects than the traditional feature extraction stitching algorithm , and the spliced image has geographic information, which has certain practical value.

Description of drawings

Figure 1 is a flow chart of the design method;

Figure 2 UAV carries POS information format;

Figure 3 is the flow chart of the overlapping area of adjacent images, in which Figure 3.1 is the velocity decomposition diagram, and Figure 3.2 is the calculation and analysis diagram of the overlapping area;

Figure 4 is a process diagram of extracting SURF feature points of adjacent images, where Figure 4.1 is a detection feature point map, Figure 4.2 is a feature point direction calculation map, and Figure 4.3 is a detected feature point matching point pair of an adjacent image;

Figure 5 is the process of correcting geographic coordinates according to SURF feature point pairs, in which Figure 5.1 shows two aerial images placed according to geographic coordinates, Figure 5.2 shows the effect of geographical coordinates placed after correction according to the method in this paper, and Figure 5.3 shows the correction of multiple pictures placement effect;

Figure 6: Effect comparison chart after fusion;

Detailed ways

The present invention will be further described below in conjunction with accompanying drawing;

The overall process of the present invention is shown in Figure 1, and is mainly divided into six major links: aerial photography image preprocessing link, aerial photography image overlapping area calculation link, image geographic coordinates calculation link, image SURF feature extraction link, and adjacent image SURF feature matching Point pair correction of geographic coordinates and image fusion. Among them, combining the adjacent image SURF feature matching points to correct the geographic coordinates of the image is the innovation layer of this method. In combination with the algorithm execution sequence, the specific implementation of the above-mentioned contents will be described below:

S1 aerial image preprocessing:

In the process of preprocessing geometric correction, the coordinate systems that need to be established in sequence are: earth coordinate system, geographic coordinate system, body coordinate system, photoelectric platform coordinate system, and camera coordinate system. The specific steps of the traditional correction are as follows: Let the transformation matrix from the photoelectric platform coordinate system to the digital camera coordinate system be R1, the transformation matrix from the body coordinate system to the photoelectric platform coordinate system is R2, and the transformation matrix from the geographic coordinate system to the body coordinate system is R3, the transformation matrix from the geodetic coordinate system to the spatial rectangular coordinate system is R4, then:

Among them, R _x (λ)R _y (λ)R _z (λ) are the rotation matrices around the x-axis, y-axis, and z-axis respectively, where α _P is the pixel scanning angle, β _P is the pixel deflection angle, χ is the azimuth, ω is the roll angle of the platform, β is the roll angle of the aircraft, θ is the heading angle, B is the longitude, L is the latitude, and H is the altitude. Therefore, the conversion from the camera coordinate system (Xc, Yc, Zc) to the earth coordinate system (Xe, Ye, Ze) can be expressed by the following relationship:

Conventional geometric correction is calculated and sampled pixel by pixel by formula (2). Considering that this project mainly stitches together a single flight belt, and the roll angle and pitch angle of the UAV platform are relatively small (when shooting approximately vertically downward), it can be considered that the image resolution is consistent, so in order to reduce the amount of calculation, improve the speed , the geometric correction only corrects the aircraft heading angle and platform azimuth angle in the horizontal direction. The rectified model looks like this:

Among them, χ is the azimuth angle, θ is the heading angle, and R ^T (χ+θ) is the transformation matrix. After the model conversion is completed, the gray value of the pixel in the new coordinate system is calculated by the bilinear difference resampling method. , so as to generate a new image matrix and complete the geometric correction.

S2: Calculation of overlapping area of aerial images

Calculating the overlapping area of adjacent images through POS information reduces the amount of calculation for splicing, and as the search area of feature points becomes smaller, the probability of mismatching may be smaller, which improves the efficiency of detection Calculate the overlapping area of adjacent images ,Specific steps are as follows:

1) Assuming that the UAV route deviates from true north at an angle of θ, and its speed is V, decompose the speed into two components in the direction of true north and true east. Denoted as V1 and V2 respectively, the velocity decomposition is shown in Figure 3.1.

2) The images of two adjacent exposure points taken continuously are numbered Pic1 and Pic2 respectively, and the interval time is t. According to the POS track file, the longitude and latitude of the center points of the two pictures are recorded as LatA, LonA, LatB, and LonB respectively, then Have:

Among them, R=6371.004 kilometers, pi is 3.1415926, C is the longitude and latitude angle between the center points of the two pictures, and L is the calculated actual distance between the two points.

3) Taking pic1 as the benchmark, the overlapping area between pic2 and pic1 is expressed as a rectangular area S as shown in Figure 3.2, and the overlapping area is expanded to define a regular area, and the vertices of the overlapping area are projected to the x and y directions to represent for:

S=(W-P1)(H-P2) (5)

Among them, W and H are the width and height of the rectangular area, and the lengths of P1 and P2 are marked in the figure.

By calculating the overlapping area, the calculation amount of splicing is reduced, and as the search area of feature points is smaller, the probability of false matching is smaller, and the detection efficiency is improved. The schematic diagram of the specific calculation overlapping area is shown in Figure 3.2. S3: Geographical coordinate calculation of the image: For the geometrically corrected aerial image, each pixel coordinate of the aerial image has a corresponding geographic coordinate, and the POS information (as shown in Figure 2) carried by the UAV during flight records the current The longitude and latitude of the aircraft in the exposure point image, as well as the altitude, roll angle, pitch angle, heading angle and other information of the aircraft can be used to calculate the geographic coordinates. The calculation process of the geographic coordinates of the image is divided into the following steps:

1) Calculate the ground resolution, the calculation formula is as follows:

Among them, GSD represents the ground resolution (m), f is the focal length of the lens (mm), P is the pixel size of the imaging sensor (mm), and H is the flying height corresponding to the UAV (m).

2) Calculate the actual ground distance of the image diagonal:

where w and h are the width and height of the image, and L is the actual distance between the diagonals of the image.

3) Calculate the geographic coordinates of the four corners of the image, and obtain the corresponding geographic coordinates of the four corners with the center point of the image as the center of the circle and a radius of L/2 according to the latitude and longitude of the center point of the image and the distance and direction angle of another point relative to the center point. The specific calculation formula is:

Where θ _i ∈ (0,2pi), Lon _a Lat _a is the latitude and longitude of the center point of the image, Ri is the equatorial radius, which is 6378137m, Rj is the polar radius, which is 6356725m, and pi is 3.1415925.

4) Geographical coordinate conversion to conversion between coordinate systems between spaces, the conversion formula is as follows:

Among them, N is the radius of curvature, and Lon, Lat, and H are the longitude, latitude, and height of any point on the image, so that the image is transformed into a space coordinate system, and the complete image coordinates are obtained so far.

S4: SURF feature extraction of image:

SURF uses an approximate Hessian matrix to detect feature points, and uses the integral image to perform convolution operations, which greatly reduces operations and thus improves feature extraction speed. The SURF descriptor consists of two main parts: detecting feature points and calculating features; the specific implementation is divided into the following steps:

1) Construct the Hessian matrix and construct the scale space.

Assuming that a point on the image is X(x, y), the matrix M at the σ scale is defined as:

Among them, L _xx is the result of the convolution of the second derivative of the Gaussian filter with X, and the meaning of L _xy is similar, and σ is the spatial scale. When the discriminant of the Hessian matrix obtains a local maximum value, it is considered to locate the position of the key point.

2) Detect feature points

In the scale space obtained above, each pixel processed by the Hessian matrix is compared with 26 points in the two-dimensional image space and the scale space neighborhood, and the key points are initially located, and then the key points with relatively weak energy are filtered out. Points and misplaced key points, and screen out the final stable feature points. The detection process is shown in Figure 4.1:

3) Determine the main direction of the feature point

With the feature point as the center and 6σ as the radius, find the Haar wavelet response in the XY direction, count the sum of the horizontal haar wavelet features and vertical haar wavelet features of all points in the 60-degree sector, and set the size of the haar wavelet to 4s. So that each sector has a corresponding value. Then rotate the 60-degree sector at a certain interval, and finally lock the direction of the maximum sector as the main direction of the feature point. A schematic diagram of the process is shown in Figure 4.2:

4) Calculate feature descriptors

Take a square frame around the feature point, and then divide the frame into 16 sub-regions, each sub-region counts the haar wavelet features of 25 pixels in the horizontal direction and vertical direction, where the horizontal and vertical directions are relative to the main direction Yes, so each feature point is a 16*4=64-dimensional vector, which will greatly speed up the matching speed during the feature matching process. An example of matching point pairs extracted from adjacent SURF feature images is shown in Figure 4.3.

S5: Use adjacent image features to match point pairs to correct geographic coordinates

Due to the low precision of POS and certain errors in geometric correction, the calculated coordinate mapping relationship has certain errors. At this time, the feature matching algorithm is used to correct the geographic coordinates of the image. The specific process is as follows:

Assuming that the geographic coordinates of image 1 are P1 (x1, y1), and the geographic coordinates of image 2 are P2 (x2, y2), after extracting image feature matching and pairing, the pixel coordinate position of the same point between two different images can be obtained. Thus, the latitude and longitude coordinates (Lon1, Lat1) (Lon2, Lat2) of the same target point in the two images are obtained, and finally, based on the geographic coordinates of the first image, the relationship between the target point in the second image and the first image is obtained. The offset between , the formula is as follows:

Then use the obtained offset to correct the geographic coordinates P2 of image 2:

Then the geographical coordinates are projected into the spatial rectangular coordinate system, so as to complete the precise registration of the image and obtain the image after corrected coordinates; Figures 5.1 and 5.2 show the effect of placing directly according to the geographical coordinates and the result of placing after correction In contrast, Figure 5.3 is the placement result of multiple images according to the corrected geographic coordinates.

S5 image fusion strategy

For the image after correcting the coordinates, since there are certain traces between each image, some strategies need to be used to solve the problem of large color transition differences between the stitching gaps, so that the stitching image is smoother and more natural, and the fade-in and fade-out are adaptive The weight fusion process is as follows:

Assuming that I1, I2, and I are image 1 before fusion, image 2 and image 3 after fusion respectively, the image fusion is completed by formula (11),

In the formula, W is the total width of the repetition of two different pictures; w is the horizontal distance between the left edge of the overlapping area and the current pixel. An example diagram of comparison before and after fusion is shown in Figure 6.

To sum up, the present invention proposes a splicing method based on the POS information of the UAV based on the characteristics of the aerial image of the UAV. First, the geographical coordinates of the four corners of the image are calculated according to the POS information, and then the overlapping area of the image is extracted. SURF features, using the positional relationship of the same feature points of adjacent images to correct the geographic coordinates, so as to complete the registration, and finally use the adaptive gradual in and out fusion algorithm to make a smooth transition of the image, and obtain a panoramic image with good visual effect , and the spliced image has geographical coordinates and has use value.

Claims (2)

1. an image mosaic method based on unmanned aerial vehicle POS information and image SURF feature combination, it is characterized in that: the method comprises the steps: S1: Aerial image preprocessing: During the aerial photography mission of the UAV, due to the influence of the aircraft's attitude, altitude, speed and the earth's rotation, the image is squeezed, twisted, stretched and offset relative to the actual position of the ground target. Carry out geometric correction to obtain remote sensing images based on the same reference projection plane; S2: Calculation of overlapping area of aerial images: Calculate the overlapping area of adjacent images according to the POS information, obtain the exposure time of the sequence image and the actual distance of the flight according to the POS information, decompose the speed in the positive direction according to the current flight angle of the aircraft, and calculate the overlapping area of adjacent captured images In order to reduce the calculation amount of splicing, and as the search area of the feature points is smaller, the probability of false matching is smaller, and the detection efficiency is improved; S3: Geographical coordinate calculation of the image: After the image is geometrically corrected, the latitude and longitude coordinates of the center point of the image are calculated according to the attitude angle of the current aircraft flight, and the ground resolution of the image is calculated according to the internal and external orientation elements of the camera, and the latitude and longitude coordinates of the four points of the image are calculated, and then the geographic coordinates are converted To the space Cartesian coordinate system, so as to project according to the space Cartesian coordinate system, to obtain the relative positional relationship between the images; S4: Extract the SURF feature of the image: The extraction of SURF features includes: constructing a Hessian matrix, generating all interest points, constructing a scale space, locating feature points, assigning the main direction of feature points, generating feature point descriptors, and matching feature points. Through these steps, the adjacent The SURF feature of the image matches the point pair; S5: Use the extracted image SURF feature point pairs to correct the geographic coordinates: Since the geographic coordinates of each image are already known, the geographic coordinates of the obtained point pairs with the same name should be the same point. If the calculated ones are different, use the first image as the reference, and correct the subsequent geographic coordinates in order to determine the position. registration to improve the accuracy of placing directly according to geographic coordinates; S6: Image fusion strategy: For the image after correcting the coordinates, since there are certain traces between each image, a strategy is used to solve the problem of color transition differences between stitching gaps, and the gradual in and gradual out adaptive weight fusion method is used to smoothly transition the gaps between images . 2. The image stitching method based on the combination of UAV POS information and image SURF features according to claim 1, characterized in that: it is mainly divided into six major links: aerial image preprocessing link, aerial image overlapping area calculation link , image geographic coordinates calculation link, image SURF feature extraction link, combined with adjacent image SURF feature matching point pair correction geographic coordinates link and image fusion link; wherein, combined with adjacent image SURF feature matching point pair correction image geographic coordinates is the method innovation layer; S1 aerial image preprocessing: In the process of preprocessing geometric correction, the coordinate systems that need to be established in turn are: earth coordinate system, geographic coordinate system, body coordinate system, photoelectric platform coordinate system, camera coordinate system; the specific steps of traditional correction are as follows: set the photoelectric platform coordinate system The transformation matrix to the digital camera coordinate system is R1, the transformation matrix from the body coordinate system to the photoelectric platform coordinate system is R2, the transformation matrix from the geographic coordinate system to the body coordinate system is R3, and the transformation matrix from the earth coordinate system to the spatial rectangular coordinate system The transformation matrix is R4, then there are: <mrow><mfenced open = "{" close = ""><mtable><mtr><mtd><mrow><mi>R</mi><mn>1</mn><mo>=</mo><mi>R</mi><mrow><mo>(</mo><msub><mi>&amp;beta;</mi><mi>P</mi></msub><mo>)</mo></mrow><mi>R</mi><mrow><mo>(</mo><msub><mi>&amp;alpha;</mi><mi>P</mi></msub><mo>)</mo></mrow></mrow></mtd></mtr><mtr><mtd><mrow><mi>R</mi><mn>2</mn><mo>=</mo><mi>R</mi><mrow><mo>(</mo><mi>&amp;chi;</mi><mo>)</mo></mrow><mi>R</mi><mrow><mo>(</mo><mi>&amp;omega;</mi><mo>)</mo></mrow></mrow></mtd></mtr><mtr><mtd><mrow><mi>R</mi><mn>3</mn><mo>=</mo><msub><mi>R</mi><mi>x</mi></msub><mrow><mo>(</mo><mi>&amp;beta;</mi><mo>)</mo></mrow><msub><mi>R</mi><mi>z</mi></msub><mrow><mo>(</mo><mi>&amp;theta;</mi><mo>)</mo></mrow><msub><mi>R</mi><mi>y</mi></msub><mrow><mo>(</mo><mn>90</mn><mo>-</mo><mi>&amp;lambda;</mi><mo>)</mo></mrow><msub><mi>R</mi><mi>x</mi></msub><mrow><mo>(</mo><mn>90</mn><mo>)</mo></mrow></mrow></mtd></mtr><mtr><mtd><mrow><mi>R</mi><mn>4</mn><mo>=</mo><msup><msub><mi>R</mi><mi>x</mi></msub><mi>T</mi></msup><mrow><mo>(</mo><mo>-</mo><mo>(</mo><mn>90</mn><mo>-</mo><mi>B</mi><mo>)</mo></mrow><mo>)</mo><msup><msub><mi>R</mi><mi>z</mi></msub><mi>T</mi></msup><mrow><mo>(</mo><mn>90</mn><mo>+</mo><mi>L</mi><mo>)</mo></mrow></mrow></mtd></mtr></mtable></mfenced><mo>-</mo><mo>-</mo><mo>-</mo><mrow><mo>(</mo><mn>1</mn><mo>)</mo></mrow></mrow> Among them, R _x (λ)R _y (λ)R _z (λ) are the rotation matrices around the x-axis, y-axis, and z-axis respectively, where α _P is the pixel scanning angle, β _P is the pixel deflection angle, χ is the azimuth, ω is the roll angle of the platform, β is the aircraft roll angle, θ is the heading angle, B is longitude, L is latitude, and H is height; therefore, from the camera coordinate system (Xc, Yc, Zc) to the earth coordinate The conversion between systems (Xe, Ye, Ze) is expressed by the following relationship: <mrow><mfenced open = "[" close = "]"><mtable><mtr><mtd><mi>X</mi><mi>e</mi></mtd></mtr><mtr><mtd><mi>Y</mi><mi>e</mi></mtd></mtr><mtr><mtd><mi>Z</mi><mi>e</mi></mtd></mtr></mtable></mfenced><mo>=</mo><msup><msub><mi>R</mi><mn>1</mn></msub><mi>T</mi></msup><msup><msub><mi>R</mi><mn>2</mn></msub><mi>T</mi></msup><msup><msub><mi>R</mi><mn>3</mn></msub><mi>T</mi></msup><msup><msub><mi>R</mi><mn>4</mn></msub><mi>T</mi></msup><mfenced open = "[" close = "]"><mtable><mtr><mtd><mi>X</mi><mi>c</mi></mtd></mtr><mtr><mtd><mi>Y</mi><mi>c</mi></mtd></mtr><mtr><mtd><mi>Z</mi><mi>c</mi></mtd></mtr></mtable></mfenced><mo>-</mo><mo>-</mo><mo>-</mo><mrow><mo>(</mo><mn>2</mn><mo>)</mo></mrow></mrow> Conventional geometric correction is calculated and sampled pixel by pixel according to formula (2); this method is aimed at splicing a single flight belt, and the roll angle and pitch angle of the UAV platform are relatively small, that is, the image resolution is considered to be high when shooting vertically downward. are consistent, the geometric correction only corrects the aircraft heading angle and platform azimuth in the horizontal direction; the correction model is as follows: <mrow><mfenced open = "[" close = "]"><mtable><mtr><mtd><mi>X</mi><mi>e</mi></mtd></mtr><mtr><mtd><mi>Y</mi><mi>e</mi></mtd></mtr><mtr><mtd><mi>Z</mi><mi>e</mi></mtd></mtr></mtable></mfenced><mo>=</mo><msup><mi>R</mi><mi>T</mi></msup><mrow><mo>(</mo><mi>&amp;chi;</mi><mo>+</mo><mi>&amp;theta;</mi><mo>)</mo></mrow><mfenced open = "[" close = "]"><mtable><mtr><mtd><mi>X</mi><mi>c</mi></mtd></mtr><mtr><mtd><mi>Y</mi><mi>c</mi></mtd></mtr><mtr><mtd><mi>Z</mi><mi>c</mi></mtd></mtr></mtable></mfenced><mo>-</mo><mo>-</mo><mo>-</mo><mrow><mo>(</mo><mn>3</mn><mo>)</mo></mrow></mrow> Among them, χ is the azimuth angle, θ is the heading angle, and R ^T (χ+θ) is the transformation matrix. After the model conversion is completed, the gray value of the pixel in the new coordinate system is calculated by the bilinear difference resampling method. , so as to generate a new image matrix and complete the geometric correction; S2: Calculation of overlapping area of aerial images 1) Assuming that the UAV flight path deviates from true north at an angle of θ, and its speed is V, decompose the speed into two components in the direction of true north and true east; denoted as V1 and V2 respectively; 2) Continuously shoot the images of two adjacent exposure points respectively numbered as Pic1 and Pic2, and the interval time is t. According to the POS track file, the longitude and latitude of the center point of the two images are recorded as LatA, LonA, LatB, and LonB respectively, then Have: <mrow><mfenced open = "{" close = ""><mtable><mtr><mtd><mrow><mi>C</mi><mo>=</mo><mi>sin</mi><mrow><mo>(</mo><mi>L</mi><mi>a</mi><mi>t</mi><mi>A</mi><mo>)</mo></mrow><mo>*</mo><mi>sin</mi><mrow><mo>(</mo><mi>L</mi><mi>a</mi><mi>t</mi><mi>B</mi><mo>)</mo></mrow><mo>*</mo><mi>cos</mi><mrow><mo>(</mo><mi>L</mi><mi>o</mi><mi>n</mi><mi>A</mi><mo>-</mo><mi>L</mo>mi><mi>o</mi><mi>n</mi><mi>B</mi><mo>)</mo></mrow><mo>+</mo><mi>cos</mi><mrow><mo>(</mo><mi>L</mi><mi>a</mi><mi>t</mi><mi>A</mi><mo>)</mo></mrow><mo>*</mo><mi>cos</mi><mrow><mo>(</mo><mi>L</mi><mi>a</mi><mi>t</mi><mi>B</mi><mo>)</mo></mrow></mrow></mtd></mtr><mtr><mtd><mrow><mi>L</mi><mo>=</mo><mi>R</mi><mo>*</mo><mi>A</mi><mi>r</mi><mi>c</mi><mi>cos</mi><mrow><mo>(</mo><mi>C</mi><mo>)</mo></mrow><mo>*</mo><mi>P</mi><mi>i</mi><mo>/</mo><mn>180</mn></mrow></mtd></mtr></mo>mtable></mfenced><mo>-</mo><mo>-</mo><mo>-</mo><mrow><mo>(</mo><mn>4</mn><mo>)</mo></mrow></mrow> Among them, R=6371.004 kilometers, pi is 3.1415926, C is the longitude and latitude angle between the center points of the two pictures, and L is the calculated actual distance between the two points; 3) Taking pic1 as the benchmark, the overlapping area between pic2 and pic1 is expressed as a rectangular area S, and the overlapping area is expanded and defined into a regular area, and the overlapping area vertices are projected to the x and y directions as follows: S=(W-P1)(H-P2) (5) Among them, W and H are the width and height of the rectangular area, and the lengths of P1 and P2 are marked in the figure; By calculating the overlapping area to reduce the amount of calculation of splicing, and as the search area of the feature points is smaller, the probability of mismatching may be smaller, improving the efficiency of detection; S3: Calculation of geographic coordinates of the image: For the geometrically corrected aerial image, each pixel coordinate of the aerial image has a corresponding geographic coordinate, and the POS information carried by the drone during flight records the latitude and longitude of the aircraft that is currently exposed to the image And the altitude, roll angle, pitch angle, and heading angle information of the aircraft are used to calculate the geographic coordinates. The calculation process of the geographic coordinates of the image is divided into the following steps: 1) Calculate the ground resolution, the calculation formula is as follows: <mrow><mi>G</mi><mi>S</mi><mi>D</mi><mo>=</mo><mfrac><mrow><mi>H</mi><mo>*</mo><mi>P</mi></mrow><mi>f</mi></mfrac><mo>-</mo><mo>-</mo><mo>-</mo><mrow><mo>(</mo><mn>6</mn><mo>)</mo></mrow></mrow> Among them, GSD represents the ground resolution, m; f is the focal length of the lens, mm; P is the pixel size of the imaging sensor, mm; H is the flight height corresponding to the UAV, m; 2) Calculate the actual ground distance of the image diagonal: <mrow><mi>L</mi><mo>=</mo><msqrt><mrow><msup><mi>w</mi><mn>2</mn></msup><mo>+</mo><msup><mi>h</mi><mn>2</mn></msup></mrow></msqrt><mo>*</mo><mi>G</mi><mi>S</mi><mi>D</mi><mo>-</mo><mo>-</mo><mo>-</mo><mrow><mo>(</mo><mn>7</mn><mo>)</mo></mrow></mrow> Where w and h are the width and height of the image, and L is the actual distance between the diagonals of the image; 3) Calculate the geographical coordinates of the four corners of the image, according to the longitude and latitude of the center point of the image and the distance and direction angle of another point relative to the center point, obtain the corresponding four corners geographical coordinates with the center point of the image as the center and a radius of L/2; the specific calculation formula is : <mrow><msub><mi>Lon</mi><mi>i</mi></msub><mo>=</mo><mfrac><mrow><mo>(</mo><mfrac><mrow><mi>L</mi><mi>sin</mi><mrow><mo>(</mo><msub><mi>&amp;theta;</mi><mi>i</mi></msub><mo>)</mo></mrow></mrow><mrow><mn>2</mn><mrow><mo>(</mo><msub><mi>R</mi><mi>j</mi></msub><mo>+</mo><mo>(</mo><msub><mi>R</mi><mi>c</mi></msub><mo>-</mo><msub><mi>R</mi><mi>j</mi></msub><mo>)</mo></mrow><mo>*</mo><mrow><mo>(</mo><mn>90</mn><mo>-</mo><msub><mi>Lat</mi><mi>a</mi></msub><mo>)</mo></mrow><mo>/</mo><mn>90</mn><mo>)</mo><mo>*</mo><mi>cos</mi><mrow><mo>(</mo><mi>L</mi><mi>a</mi><mi>t</mi><mi>a</mi><mo>)</mo></mrow></mrow></mfrac><mo>+</mo><msub><mi>Lon</mi><mi>a</mi></msub><mo>)</mo><mo>*</mo><mn>180</mn></mrow><mrow><mi>p</mi><mi>i</mi></mrow></mfrac><mo>-</mo><mo>-</mo><mo>-</mo><mrow><mo>(</mo><mn>8</mn><mo>)</mo></mrow></mrow> <mrow><msub><mi>Lat</mi><mi>i</mi></msub><mo>=</mo><mfrac><mrow><mo>(</mo><mfrac><mrow><mi>L</mi><mi></mi><mi>c</mi><mi>o</mi><mi>s</mi><mrow><mo>(</mo><msub><mi>&amp;theta;</mi><mi>i</mi></msub><mo>)</mo></mrow></mrow><mrow><mn>2</mn><mrow><mo>(</mo><msub><mi>R</mi><mi>j</mi></msub><mo>+</mo><mo>(</mo><msub><mi>R</mi><mi>c</mi></msub><mo>-</mo><msub><mi>R</mi><mi>j</mi></msub><mo>)</mo></mrow><mo>*</mo><mrow><mo>(</mo><mn>90</mn><mo>-</mo><msub><mi>Lat</mi><mi>a</mi></msub><mo>)</mo></mrow><mo>/</mo><mn>90</mn><mo>)</mo></mrow></mfrac><mo>+</mo><msub><mi>Lat</mi><mi>a</mi></msub><mo>)</mo><mo>*</mo><mn>180</mn></mrow><mrow><mi>p</mi><mi>i</mi></mrow></mfrac><mo>-</mo><mo>-</mo><mo>-</mo><mrow><mo>(</mo><mn>9</mn><mo>)</mo></mrow></mrow> Where θ _i ∈ (0,2pi), Lon _a Lat _a is the latitude and longitude of the center point of the image, Ri is the equatorial radius, which is 6378137m, Rj is the polar radius, which is 6356725m, and pi is 3.1415925; 4) Geographical coordinate conversion to conversion between coordinate systems between spaces, the conversion formula is as follows: <mrow><mfenced open = "[" close = "]"><mtable><mtr><mtd><msub><mi>X</mi><mi>s</mi></msub></mrow>mtd></mtr><mtr><mtd><msub><mi>Y</mi><mi>s</mi></msub></mtd></mtr><mtr><mtd><msub><mi>Z</mi><mi>s</mi></msub></mtd></mtr></mtable></mfenced><mo>=</mo><mfenced open = "[" close = "]"><mtable><mtr><mtd><mo>(</mo><mi>N</mi><mo>+</mo><mi>H</mi><mo>)</mo><mi>c</mi><mi>o</mi><mi>s</mi><mo>(</mo><mi>L</mi><mi>o</mi><mi>n</mi><mo>)</mo><mi>c</mi><mi>o</mi><mi>s</mi><mo>(</mo><mi>L</mi><mi>a</mi><mi>t</mi><mo>)</mo></mtd></mtr><mtr><mtd><mo>(</mo><mi>N</mi><mo>+</mo><mi>H</mi><mo>)</mo><mi>c</mi><mi>o</mi><mi>s</mi><mo>(</mo><mi>L</mi><mi>o</mi><mi>n</mi><mo>)</mo><mi>sin</mi><mo>(</mo><mi>L</mi><mi>a</mi><mi>t</mi><mo>)</mo>mo></mtd></mtr><mtr><mtd><mo>&amp;lsqb;</mo><mi>N</mi><mo>(</mo><mn>1</mn><mo>-</mo><msup><mi>e</mi><mn>2</mn></msup><mo>)</mo><mo>+</mo><mi>H</mi><mo>&amp;rsqb;</mo><mi>sin</mi><mi></mi><mi>L</mi><mi>o</mi><mi>n</mi></mtd></mtr></mtable></mfenced><mo>-</mo><mo>-</mo><mo>-</mo><mrow><mo>(</mo><mn>10</mn><mo>)</mo></mrow></mrow> Among them, N is the radius of curvature, Lon, Lat, and H are the longitude, latitude and height of any point on the image, and the image is converted to the space coordinate system, so far the complete image coordinates are obtained; S4: SURF feature extraction of image: SURF uses an approximate Hessian matrix to detect feature points, and uses integral images for convolution operations, reducing operations and improving feature extraction speed; SURF descriptors include two parts: detecting feature points and calculating features; the specific implementation is divided into the following steps: 1) Construct the Hessian matrix and construct the scale space; Assuming that a point on the image is X(x, y), the matrix M at the σ scale is defined as: <mrow><mi>M</mi><mrow><mo>(</mo><mi>x</mi><mo>,</mo><mi>&amp;sigma;</mi><mo>)</mo></mrow><mo>=</mo><mfenced open = "[" close = "]"><mtable><mtr><mtd><mrow><msub><mi>L</mi><mrow><mi>x</mi><mi>x</mi></mrow></msub><mrow><mo>(</mo><mi>x</mi><mo>,</mo><mi>&amp;sigma;</mi><mo>)</mo></mrow></mrow></mtd><mtd><mrow><msub><mi>L</mi><mrow><mi>x</mi><mi>y</mi></mrow></msub><mrow><mo>(</mo><mi>x</mi><mo>,</mo><mi>&amp;sigma;</mi><mo>)</mo></mrow></mrow></mtd></mtr><mtr><mtd><mrow><msub><mi>L</mi><mrow><mi>x</mi><mi>y</mi></mrow></msub><mrow><mo>(</mo><mi>x</mi><mo>,</mo><mi>&amp;sigma;</mi><mo>)</mo></mrow></mrow></mtd><mtd><mrow><msub><mi>L</mi><mrow><mi>y</mi><mi>y</mi></mrow></msub><mrow><mo>(</mo><mi>x</mi><mo>,</mo><mi>&amp;sigma;</mi><mo>)</mo></mrow></mrow></mtd></mtr></mtable></mfenced><mo>-</mo><mo>-</mo><mo>-</mo><mrow><mo>(</mo><mn>11</mn><mo>)</mo></mrow></mrow> Among them, _{Lxx is the result of Gaussian filter second-order derivative and X convolution, Lxy ,} etc. have similar meanings, and σ is the spatial scale; when the discriminant of the Hessian matrix obtains a local maximum value, it is considered to be located at the position of the key point; 2) Detect feature points In the obtained scale space, each pixel processed by the Hessian matrix is compared with 26 points in the two-dimensional image space and the scale space neighborhood, and the key points are initially located, and then the key points with weak energy are removed by filtering and The key points of wrong positioning are screened out to the final stable feature points; 3) Determine the main direction of the feature point With the feature point as the center and 6σ as the radius, find the Haar wavelet response in the XY direction, count the sum of the horizontal haar wavelet features and vertical haar wavelet features of all points in the 60-degree sector, and set the size of the haar wavelet to 4s. So that each sector has a corresponding value; then rotate the 60-degree sector at a certain interval, and finally lock the direction of the maximum sector as the main direction of the feature point; 4) Calculate feature descriptors Take a square frame around the feature point, and then divide the frame into 16 sub-regions, each sub-region counts the haar wavelet features of 25 pixels in the horizontal direction and vertical direction, where the horizontal and vertical directions are relative to the main direction , so each feature point is a 16*4=64-dimensional vector; S5: Use adjacent image features to match point pairs to correct geographic coordinates Due to the low precision of POS and certain errors in geometric correction, the calculated coordinate mapping relationship has certain errors. At this time, the feature matching algorithm is used to correct the geographic coordinates of the image. The specific process is as follows: Assuming that the geographic coordinates of image 1 are P1 (x1, y1), and the geographic coordinates of image 2 are P2 (x2, y2), after extracting image feature matching and pairing, the pixel coordinate position of the same point between two different images can be obtained. Thus, the latitude and longitude coordinates (Lon1, Lat1) (Lon2, Lat2) of the same target point in the two images are obtained, and finally, based on the geographic coordinates of the first image, the relationship between the target point in the second image and the first image is obtained. The offset between , the formula is as follows: <mrow><mfenced open = "{" close = ""><mtable><mtr><mtd><mrow><mi>x</mi><mo>=</mo><mi>L</mi><mi>o</mi><mi>n</mi><mn>2</mn><mo>-</mo><mi>L</mi><mi>o</mi><mi>n</mi><mn>1</mn></mrow></mtd></mtr><mtr><mtd><mrow><mi>y</mi><mo>=</mo><mi>L</mi><mi>a</mi><mi>t</mi><mn>2</mn><mo>-</mo><mi>L</mi><mi>a</mi><mi>t</mi><mn>1</mn></mrow></mtd></mtr></mtable></mfenced><mo>-</mo><mo>-</mo><mo>-</mo><mrow><mo>(</mo><mn>12</mn><mo>)</mo></mrow></mrow> Then use the obtained offset to correct the geographic coordinates P2 of image 2: <mrow><mfenced open = "{" close = ""><mtable><mtr><mtd><mrow><mi>x</mi><mn>2</mn><mo>=</mo><mi>x</mi><mn>2</mn><mo>+</mo><mi>x</mi></mrow></mtd></mtr><mtr><mtd><mrow><mi>y</mi><mn>2</mn><mo>=</mo><mi>y</mi><mn>2</mn><mo>+</mo><mi>y</mi></mrow></mtd></mtr></mtable></mfenced><mo>-</mo><mo>-</mo><mo>-</mo><mrow><mo>(</mo><mn>13</mn><mo>)</mo></mrow></mrow> Then the geographic coordinates are projected into the space Cartesian coordinate system, so as to complete the precise registration of the image and obtain the image after correcting the coordinates; S5 image fusion strategy For the image after correcting the coordinates, since there are certain traces between each image, some strategies need to be used to solve the problem of large color transition differences between the stitching gaps, so that the stitching image is smoother and more natural, and the fade-in and fade-out are adaptive The weight fusion process is as follows: Assuming that I1, I2, and I are image 1 before fusion, image 2 and image 3 after fusion respectively, the image fusion is completed by formula (11), In the formula, W is the total width of the repetition of two different pictures; w is the horizontal distance between the left edge of the overlapping area and the current pixel.