An unmanned aerial vehicle image matching method based on spatio-temporal change and neighborhood category outlier

Abstract

The application discloses a UAV image matching method based on space-time change and neighborhood category outlier, and belongs to the technical field of UAV image processing and registration. The application solves the problem of low image registration accuracy caused by inaccurate estimation of homographic transformation matrix in the prior art. The application uses a space-time change description vector construction module to extract space-time change description vectors from the flight azimuth difference, the pitch angle difference, the spatial displacement and the shooting time interval of two adjacent images. For each matching point pair, neighborhood space-time change category vectors are constructed according to the displacement direction deviation, the displacement amplitude ratio of the neighboring matching point pairs and the alignment degree with the global space-time change description. Then, the neighborhood space-time change category vectors of all matching point pairs are subjected to K-means clustering. According to the clustering result, the error matching point pairs are determined and removed. Finally, the homographic transformation matrix is accurately estimated according to the remaining correct matching point pair set, and the registered image is obtained. The application can be applied to UAV image matching.

Description

Technical Field

[0001] This invention belongs to the field of UAV image processing and registration, specifically relating to a UAV image matching method based on spatiotemporal changes and neighborhood category outliers. Background Technology

[0002] UAV multi-angle image stitching and 3D reconstruction are fundamental technologies in precision agriculture, ecological surveys, and topographic mapping. The core process involves extracting feature points between adjacent images and performing initial matching. Based on the matched point pairs, the geometric transformation relationship between the images is calculated, leading to image registration and stitching. When the subject is a man-made target such as a building or road, the images contain rich edges and regular textures, resulting in high feature matching accuracy. The homography transformation matrix estimated from the matched point set is reliable, and the image registration result is accurate. However, in agricultural planting areas, garden landscapes, and natural vegetation communities, the plant canopy texture is complex and highly repetitive, with very few distinct edge features. Feature matching algorithms generate a large number of incorrect matching points. Once these incorrect matching points participate in geometric transformation estimation, they directly cause serious deviations in the homography transformation matrix estimation results, image registration failure, stitching misalignment, and 3D point cloud deformation, severely limiting the application effectiveness of UAV remote sensing data in vegetation monitoring.

[0003] Existing methods primarily perform matching and filtering at the visual feature level, which has limited effectiveness in scenarios with extremely high mismatch rates, such as those involving plant canopies. Current UAV image registration mainly relies on the following technical means:

[0004] I. A matching-filtered registration method based on the distance ratio of feature descriptors uses the ratio of the nearest neighbor distance to the second nearest neighbor distance of feature descriptors such as SIFT or ORB. Matches with a ratio less than a threshold (e.g., 0.75) are retained, while those higher are discarded. The homography transformation matrix is ​​then estimated based on the retained matching point pairs to complete image registration. The main limitations of this method are: the distance ratio threshold is a globally fixed value. In areas with highly repetitive textures, such as plant canopies, many incorrectly matched descriptors are very close to the correct matches. Therefore, the distance ratio cannot effectively distinguish incorrect matching point pairs. Furthermore, this method judges each point independently, without utilizing the spatial consistency constraint between adjacent matching points. The remaining incorrectly matched point pairs after filtering still lead to significant deviations in the homography transformation matrix estimation, ultimately resulting in obvious misalignment or even failure in image registration.

[0005] II. A global geometric model-based registration method using random sampling consistency employs the RANSAC algorithm to randomly sample matching point pairs from the set of matching point pairs. It then estimates the homography transformation matrix or fundamental matrix based on the sampled matching point pairs. Matching point pairs with reprojection errors exceeding a threshold are identified as outliers and removed. Finally, the inlier set is used to estimate the final transformation and complete the registration. The main limitations of this method are: the global single geometric model assumes the scene is planar or near-planar. In vegetation scenes with significant differences in canopy height, the single model cannot accurately describe the spatial relationships of all correctly matched points, leading to some correctly matched points being misclassified as outliers. When the mismatch ratio exceeds 50%, the RANSAC algorithm has difficulty converging, and the mismatch ratio in vegetation scenes often reaches or exceeds this level, ultimately resulting in inaccurate transformation estimation and low-quality registered images.

[0006] III. A registration method based on visual consistency matching verification involves extracting local image patches around the matching point for cross-correlation or template matching verification. Matches with consistency below a threshold are discarded, and the verified matching points drive image registration. The main limitations of this method are: similar texture patches in the plant canopy appear repeatedly in different locations, meaning highly visually consistent matches are not necessarily spatially correct; this method relies entirely on image visual information, failing to utilize external parameters such as the drone's shooting direction and flight time, ignoring the spatiotemporal constraints between consecutively captured images, resulting in a large number of erroneous matches remaining and participating in transform estimation, severely degrading the final registration accuracy.

[0007] IV. A registration method based on statistical threshold-based displacement filtering calculates the displacement vectors of all matching points from image A to image B. Using the mean and standard deviation of the displacement as a benchmark, matching point pairs whose displacement deviates from the mean by more than a certain multiple of the standard deviation are discarded. The homography transformation matrix is ​​estimated based on the remaining matching point pairs after filtering, and registration is completed. The main limitations of this method are: the global statistical mean is not uniform when the image is rotated or the viewpoint changes; the displacement of matching points in the edge region is different from that in the center region; and a fixed multiple of the standard deviation threshold cannot adapt to the spatial distribution differences of displacement. Furthermore, it does not consider the local neighborhood consistency of matching point pairs, and cannot distinguish between isolated random errors and locally clustered systematic errors. Residual clustered errors can still cause local deformation in the registered image.

[0008] In summary, existing UAV image registration methods suffer from several key shortcomings, including the inability of global thresholds to adapt to local texture differences, the inapplicability of single geometric models to multi-layered canopy scenes, the neglect of flight spatiotemporal parameters by purely visual methods, and the failure of global statistics to consider neighborhood spatial consistency. These shortcomings result in a large number of erroneous matching point pairs remaining in the cleaned matching point pair set, leading to inaccurate estimation of the homography transformation matrix and low image registration accuracy. Therefore, it is urgent to propose a new UAV image matching method to solve the above problems. Summary of the Invention

[0009] The purpose of this invention is to solve the problem of low image registration accuracy caused by inaccurate estimation of homography transformation matrix in existing methods, and to propose a UAV image matching method based on spatiotemporal changes and neighborhood category outliers.

[0010] The technical solution adopted by this invention to solve the above-mentioned technical problems is: a UAV image matching method based on spatiotemporal changes and neighborhood category outliers, the method specifically includes the following steps:

[0011] Step 1: Acquire images A and B to be registered, and record the flight parameters of the UAV when capturing images A and B respectively. The flight parameters include the azimuth angle, pitch angle, flight altitude, X coordinate of the UAV center, Y coordinate of the UAV center, and capture time when capturing the images.

[0012] Perform initial matching on image A and image B to obtain an initial matching point pair list MatchList. Each record in the MatchList includes the following fields:

[0013] The MPID field is the number of the matching point pair;

[0014] The MPxA field represents the pixel coordinates of the matched point pair in the width direction within image A.

[0015] The MPyA field represents the pixel coordinates of the matching point pair in image A in the height direction.

[0016] The MPxB field represents the pixel coordinates of the matched point pair in the width direction within image B.

[0017] The MPyB field represents the pixel coordinates of the matching point pair in the height direction within image B.

[0018] The MPdx field represents the value of MPxB - MPxA;

[0019] The MPdy field represents the value of MPyB-MPyA;

[0020] And record the pixel data matrices corresponding to image A and image B respectively;

[0021] Step 2: Based on the initial matching point pair list MatchList, the pixel data matrices corresponding to image A and image B, and using the spatiotemporal change description vector construction module SKBQModel, the neighborhood category outlier matching point detection module LYFCModel, and the clean point set homography transformation registration module JHPYModel, the registered image is obtained.

[0022] The beneficial effects of this invention are:

[0023] This invention proposes a spatiotemporal change description vector construction module, a neighborhood category outlier matching point detection module, and a cleaned point set homography transformation registration module. The spatiotemporal change description vector construction module extracts an 8-dimensional spatiotemporal change description vector from the flight azimuth difference, pitch difference, spatial displacement, and shooting time interval of two adjacent images to characterize the global attitude change features between images. For each matching point pair, a 12-dimensional neighborhood spatiotemporal change category vector is constructed based on the displacement direction deviation, displacement amplitude ratio, and alignment degree with the global spatiotemporal change description of its nearest matching point pairs. Then, K-means clustering is performed on the neighborhood spatiotemporal change category vectors of all matching point pairs. Matching point pairs that do not belong to the dominant cluster and whose distance to the center of the dominant cluster exceeds a threshold are judged as incorrect matching point pairs and removed. The cleaned point set homography transformation registration module accurately estimates the homography transformation matrix based on the remaining correct matching point pairs and performs perspective transformation resampling on image B. Finally, a registered image that is accurately aligned with image A is obtained, realizing automatic matching and registration of UAV images for complex texture scenes such as plants. Attached Figure Description

[0024] Figure 1 This is a flowchart of a UAV image matching method based on spatiotemporal changes and neighborhood category outliers according to the present invention. Detailed Implementation

[0025] Specific implementation method one: Combining Figure 1 This embodiment describes a UAV image matching method based on spatiotemporal changes and neighborhood category outliers. The method specifically includes the following steps:

[0026] Step 1: Acquire images A and B to be registered, and record the flight parameters of the UAV when capturing images A and B respectively. The flight parameters include the azimuth angle, pitch angle, flight altitude, X coordinate of the UAV center, Y coordinate of the UAV center, and capture time when capturing the images.

[0027] Using existing methods, initial matching is performed on images A and B to obtain an initial matching point pair list MatchList. Each record in the MatchList includes the following fields:

[0028] The MPID field is the number of the matching point pair (the number is an integer, starting from 1).

[0029] The MPxA field represents the pixel coordinates of the matched point pair in the width direction within image A.

[0030] The MPyA field represents the pixel coordinates of the matching point pair in image A in the height direction.

[0031] The MPxB field represents the pixel coordinates of the matched point pair in the width direction within image B.

[0032] The MPyB field represents the pixel coordinates of the matching point pair in the height direction within image B.

[0033] The MPdx field represents the value of MPxB - MPxA;

[0034] The MPdy field represents the value of MPyB-MPyA;

[0035] The pixel data matrices corresponding to image A and image B are recorded respectively. The number of rows in the pixel data matrix is ​​the height ImgWidth of the image, and the number of columns in the pixel data matrix is ​​the width ImgWidth of the image. Each element value in the matrix is ​​the pixel value at the corresponding position in the image. The default value of ImgWidth is 4000 and the default value of ImgHeight is 3000.

[0036] Step 2: Based on the initial matching point pair list MatchList, the pixel data matrices corresponding to image A and image B, and using the spatiotemporal change description vector construction module SKBQModel, the neighborhood category outlier matching point detection module LYFCModel, and the clean point set homography transformation registration module JHPYModel, the registered image is obtained.

[0037] It should be noted that the units for azimuth and pitch angles are degrees, with a range of 0 to 360; the units for flight altitude and UAV center coordinates are meters, and the UAV center coordinates refer to the coordinates of the UAV in a Cartesian coordinate system; the unit for shooting time is seconds, starting from the start of the mission.

[0038] This invention has the following outstanding features:

[0039] On the one hand, the spatiotemporal change description vector construction module does not rely on image visual features for matching quality judgment. Instead, it starts from external spatiotemporal parameters such as the azimuth difference, pitch difference, and shooting time interval of the UAV flight, and constructs a global attitude change description between images in the form of triangular coding and exponential decay function. This provides a spatial consistency reference benchmark for the neighborhood category outlier matching point detection module that is not affected by texture complexity. On the other hand, the neighborhood category outlier matching point detection module does not use a single global threshold or a single geometric model to judge the matching points independently. Instead, it constructs a 12-dimensional category vector for each matching point pair, which includes neighborhood displacement deviation, amplitude ratio, and global alignment. It uses K-means clustering to allow all matching point pairs to spontaneously form behavioral categories. The dominant category represents the correct matching group. It does not rely on manually set absolute thresholds. This significantly reduces the false rejection of correct matches while improving the false match detection rate, so that the matching point set entering the subsequent registration stage is fully purified.

[0040] On the other hand, the homography transformation registration module takes a high-quality set of correctly matched point pairs that has been purified by both spatiotemporal change constraints and neighborhood category outliers as input. It constructs the coefficient matrix using the direct linear transformation least squares method and estimates the homography transformation matrix through singular value decomposition. Then, it resamples the image B by pixel-by-pixel bilinear interpolation using perspective inverse transformation to achieve image registration. Since erroneous point pairs in the input matching point pair set have been fully eliminated, the estimation accuracy of the homography transformation matrix is ​​significantly improved. Finally, the alignment accuracy and reliability of the registered image are greatly improved.

[0041] With these two features, the present invention can effectively complete the automatic matching and registration of UAV images in complex texture scenes such as plants, significantly improving the accuracy and reliability of image stitching and 3D reconstruction.

[0042] Specific Implementation Method Two: This implementation method is a further limitation of Specific Implementation Method One. The spatiotemporal change description vector construction module SKBQModel is used to extract global attitude change features from the azimuth difference, pitch difference, spatial displacement and shooting time interval of the two images.

[0043] The input to the spatiotemporal change description vector construction module SKBQModel is the flight parameters taken by the UAV when it captures images A and B, and the output is the spatiotemporal change description vector SKBQVec.

[0044] The other steps and parameters are the same as in Specific Implementation Method 1.

[0045] Specific Implementation Method 3: This implementation method is a further limitation of Specific Implementation Method 2. The working process of the spatiotemporal transition description vector construction module SKBQModel is as follows:

[0046] S101. Initialize the spatiotemporal transition description vector SKBQVec = an 8-dimensional all-zero vector;

[0047] S102, Calculate the azimuth difference in radians SKBQTemp1:

[0048] SKBQTemp1=(ImgB_Azi-ImgA_Azi)×π / 180

[0049] Where π is the mathematical constant pi, ImgA_Azi is the azimuth angle when the drone captures image A, and ImgB_Azi is the azimuth angle when the drone captures image B.

[0050] Calculate the pitch angle difference in radians using SKBQTemp2:

[0051] SKBQTemp2=(ImgB_Ele-ImgA_Ele)×π / 180

[0052] Where ImgA_Ele is the pitch angle when the drone captures image A, and ImgB_Ele is the pitch angle when the drone captures image B;

[0053] S103, let SKBQVec[1]=sin(SKBQTemp1);

[0054] Wherein, SKBQVec[1] is the first element in the spatiotemporal change description vector SKBQVec, and sin is the sine function;

[0055] Let SKBQVec[2] = cos(SKBQTemp1);

[0056] Wherein, cos is the cosine function, and SKBQVec[2] is the second element in the spatiotemporal change description vector SKBQVec;

[0057] Let SKBQVec[3]=sin(SKBQTemp2), SKBQVec[4]=cos(SKBQTemp2);

[0058] Among them, SKBQVec[3] is the third element in the spatiotemporal change description vector SKBQVec, and SKBQVec[4] is the fourth element in the spatiotemporal change description vector SKBQVec.

[0059] S104. Calculate the normalized horizontal displacement SKBQTemp3 and the normalized vertical displacement SKBQTemp4:

[0060]

[0061] Where ImgA_CX is the X coordinate of the drone's center when image A is captured, ImgB_CX is the X coordinate of the drone's center when image B is captured, ABS is the absolute value, ImgA_CY is the Y coordinate of the drone's center when image A is captured, and ImgB_CY is the Y coordinate of the drone's center when image B is captured.

[0062]

[0063] Let SKBQVec'[5]=SKBQTemp3, SKBQVec[6]=SKBQTemp4;

[0064] Among them, SKBQVec'[5] is the 5th element in the spatiotemporal change description vector SKBQVec, and SKBQVec[6] is the 6th element in the spatiotemporal change description vector SKBQVec;

[0065] S105, Calculate the normalized value of height change SKBQTemp5:

[0066] SKBQTemp5=(ImgB_Alt-ImgA_Alt) / (ABS(ImgB_Alt-ImgA_Alt)+5.0)

[0067] Where ImgA_Alt is the altitude when the drone captured image A, and ImgB_Alt is the altitude when the drone captured image B;

[0068] Let SKBQVec[7] = SKBQTemp5;

[0069] Among them, SKBQVec[7] is the 7th element in the spatiotemporal change description vector SKBQVec;

[0070] S106, Calculate the time interval attenuation SKBQTemp6:

[0071] SKBQTemp6=ABS(ImgB_Time-ImgA_Time)

[0072] Where ImgA_Time is the time when image A was captured, and ImgB_Time is the time when image B was captured;

[0073] Let SKBQVec[8]=exp(-SKBQTemp6 / TimeTau);

[0074] Where exp represents an exponential function with the natural constant e as the base, TimeTau is the time decay constant in seconds, and SKBQVec[8] is the 8th element in the spatiotemporal change description vector SKBQVec;

[0075] S107, Calculate the azimuth-elevation coupling cross term SKBQTemp7:

[0076] SKBQTemp7=sin(SKBQTemp1)×cos(SKBQTemp2)×0.5

[0077] The 5th element in the spacetime transition description vector SKBQVec is updated as follows:

[0078] SKBQVec[5]=SKBQVec'[5]+SKBQTemp7×0.3

[0079] Wherein, SKBQVec[5] represents the 5th element in the final spatiotemporal transition description vector SKBQVec;

[0080] S108. Output the spatiotemporal change description vector SKBQVec as the result of the SKBQModel module.

[0081] The other steps and parameters are the same as in Specific Implementation Method Two.

[0082] Specific Implementation Method Four: This implementation method is a further limitation of Specific Implementation Method Three. The neighborhood category outlier matching point detection module LYFCModel is used to construct a neighborhood spatiotemporal change category vector for each matching point pair, and to discover and mark erroneous matching point pairs that do not belong to the dominant category in a clustering manner.

[0083] The neighborhood category outlier matching point detection module LYFCModel takes the initial matching point pair list MatchList and the spatiotemporal transition description vector SKBQVec as input, and outputs the purified valid matching point pair list JHPYValidList.

[0084] The other steps and parameters are the same as in Specific Implementation Method 3.

[0085] Specific Implementation Method Five: This implementation method further defines Specific Implementation Method Four. The working process of the neighborhood category outlier matching point detection module LYFCModel is as follows:

[0086] S201. Initialize the error point pair list LYFCErrorList = empty list;

[0087] S202. Calculate the global mean of the displacement vector in the width direction and the global mean of the displacement vector in the height direction, and then calculate the global displacement mean amplitude LYFCMeanMag based on the global mean of the displacement vector in the width direction and the global mean of the displacement vector in the height direction.

[0088] S203. Establish a point-by-point category vector calculation counter LYFCCounter1=1;

[0089] S204. Retrieve the LYFCCounter1 record from the MatchList and obtain the MPxA, MPyA, MPdx, and MPdy fields of the LYFCCounter1 record.

[0090] Calculate the displacement magnitude LYFCCurMag of the matching point pair corresponding to the LYFCCounter1 record:

[0091]

[0092] Calculate the displacement direction LYFCCurAngle of the matching point pair corresponding to the LYFCCounter1 record:

[0093]

[0094] Where arctan2 is the arctangent function for the four quadrants (returns radians);

[0095] S205. Calculate the Euclidean distance between each other matching point in image A and the matching point corresponding to the LYFCCounter1 record in image A, and store the NbK matching points with the smallest distances into the neighborhood list LYFCNbList.

[0096] Where NbK is the number of nearest neighbors, and the default value of NbK is 6;

[0097] S206. Calculate the mean displacement of the neighborhood in the width direction (LYFCNbMeanDx) and the mean displacement in the height direction (LYFCNbMeanDy) of the corresponding matching point of the LYFCCounter1 record in image A, and then calculate the magnitude of the mean displacement of the neighborhood (LYFCNbMeanMag) of the corresponding matching point of the LYFCCounter1 record in image A.

[0098] And calculate the mean of the neighborhood displacement direction LYFCNbMeanAngle and the variance of the neighborhood displacement direction LYFCNbAngVar of the corresponding matching point in image A for the LYFCCounter1 record;

[0099] S207. Construct the 12-dimensional neighborhood category vector LYFCVec for the matching point pair corresponding to the LYFCCounter1 record in the list MatchList;

[0100] S208. Assign the neighborhood category vector LYFCVec to the MPNbVec field of the LYFCCounter1 record in the list MatchList.

[0101] Let LYFCCounter1 = LYFCCounter1 + 1;

[0102] If LYFCCounter1 is less than or equal to MPTotal, proceed to S204; otherwise, proceed to S209.

[0103] S209. Set the number of clusters to ClusterNum, perform K-means clustering (K-Means algorithm, with a maximum number of iterations of 50) on the MPNbVec field of all records in the list MatchList, obtain the cluster label MPCluster for each matching pair based on the clustering results, and obtain the matrix LYFCCenters composed of the cluster centers of each cluster (the matrix has the dimension of ClusterNum rows × 12 columns).

[0104] S210. Count the number of matching point pairs contained in each cluster, select the cluster with the most matching point pairs as the dominant cluster DominantCluster, and denote the cluster center vector corresponding to the dominant cluster as DominantCenter.

[0105] S211. Establish an outlier detection counter LYFCCounter3=1;

[0106] S212. Retrieve the LYFCCounter3rd record from the MatchList list;

[0107] If the clustering label MPCluster of the matching point pair corresponding to the LYFCCounter3 record is equal to DominantCluster, then proceed to S213;

[0108] If the cluster label MPCluster of the matching point pair corresponding to the LYFCCounter3 record is not equal to DominantCluster, then calculate the Euclidean distance LYFCDistToDom between the MPNbVec field of the matching point pair corresponding to the LYFCCounter3 record and DominantCenter (that is, subtract the corresponding values ​​of each dimension, then calculate the sum of squares of each difference, and then calculate the square root of the sum).

[0109] The distance from the MPNbVec field to the DominantCenter for each matching point pair within the dominant cluster is calculated, and the standard deviation LYFCDomStd of the calculated distance is then calculated (cached and reused after the first calculation).

[0110] If LYFCDistToDom is greater than LYFCDomStd×OutlierMult, then set the error flag MPIsError=1 for the matching point pair corresponding to the LYFCCounter3 record; add the matching point pair corresponding to the LYFCCounter3 record to the error point pair list LYFCErrorList;

[0111] If LYFCDistToDom is less than or equal to LYFCDomStd×OutlierMult, then set the error flag MPIsError=0 for the matching point pair corresponding to the LYFCCounter3 record;

[0112] Where OutlierMult is the outlier distance factor, and the default value of OutlierMult is 2.0;

[0113] S213. Set the error flag MPIsError=0 for the matching point pair corresponding to the LYFCCounter3 record, and let LYFCCounter3=LYFCCounter3+1;

[0114] If LYFCCounter3 is less than or equal to MPTotal, go to S212; otherwise, go to S214.

[0115] S214, S2 ends.

[0116] The other steps and parameters are the same as in Specific Implementation Method Four.

[0117] Specific Implementation Method Six: This implementation method is a further limitation of Specific Implementation Method Five. The specific process of S202 is as follows:

[0118] Step 1: Establish a global mean temporary storage value for the width direction displacement vector LYFCMeanDx'=0 and a global mean temporary storage value for the height direction displacement vector LYFCMeanDy'=0;

[0119] And establish a summation counter LYFCCounter0=1;

[0120] Step 2: Set LYFCMeanDx' = LYFCMeanDx' + the value of the MPdx field of the LYFCCounter0th record in the list MatchList;

[0121] Let LYFCMeanDy' = LYFCMeanDy' + the value of the MPdy field of the LYFCCounter0th record in the list MatchList;

[0122] Step 3: Set LYFCCounter0 = LYFCCounter0 + 1;

[0123] If LYFCCounter0 is less than or equal to the total number of records in the MatchList (MPTotal), proceed to step 2; otherwise, proceed to step 4.

[0124] Step 4: Final global mean of the width direction displacement vector LYFCMeanDx:

[0125] LYFCMeanDx=LYFCMeanDx' / MPTotal

[0126] The final global mean of the height displacement vector is LYFCMeanDy.

[0127] LYFCMeanDy=LYFCMeanDy' / MPTotal

[0128] Step 5: Calculate the global mean displacement magnitude LYFCMeanMag based on LYFCMeanDx and LYFCMeanDy:

[0129]

[0130] The other steps and parameters are the same as in Specific Implementation Method 5.

[0131] Specific Implementation Method Seven: This implementation method is a further limitation of Specific Implementation Method Six. The specific process of S206 is as follows:

[0132] Step 1) Initialize the temporary storage of the average displacement of the neighborhood in the width direction LYFCNbMeanDx'=0, the temporary storage of the average displacement of the neighborhood in the height direction LYFCNbMeanDy'=0, and the neighborhood displacement direction list LYFCNbAngles=empty list;

[0133] Set the neighborhood traversal counter LYFCCounter2=1;

[0134] Step 2) Extract the MPdx and MPdy fields corresponding to the LYFCCounter2 neighbor point in the list LYFCNbList in the list MatchList;

[0135] Let LYFCNbMeanDx'=LYFCNbMeanDx'+ ;

[0136] LYFCNbMeanDy'=LYFCNbMeanDy'+ ;

[0137] Calculate the displacement direction of the LYFCCounter2th neighboring point in the list LYFCNbList. ,Will Add to list LYFCNbAngles;

[0138]

[0139] in, The value of the MPdx field for the LYFCCounter2 neighboring point. The value of the MPdy field for the LYFCCounter2 neighboring point. The displacement direction of the LYFCCounter2th neighboring point;

[0140] Step 3) Set LYFCCounter2 = LYFCCounter2 + 1;

[0141] If LYFCCounter2 is less than or equal to NbK, proceed to step 2; otherwise, proceed to step 4.

[0142] Step 4) Calculate the mean displacement of the neighborhood in the width direction, LYFCNbMeanDx:

[0143] LYFCNbMeanDx=LYFCNbMeanDx' / NbK

[0144] Calculate the mean displacement of the neighborhood in the height direction, LYFCNbMeanDy:

[0145] LYFCNbMeanDy=LYFCNbMeanDy' / NbK

[0146] Calculate the mean displacement magnitude of the neighborhood LYFCNbMeanMag:

[0147]

[0148] Calculate the mean displacement direction of the neighborhood LYFCNbMeanAngle:

[0149]

[0150] Step 5) Calculate the neighborhood displacement direction variance LYFCNbAngVar:

[0151]

[0152] in, This refers to the first element in the list LYFCNbAngles. Values.

[0153] The other steps and parameters are the same as in Specific Implementation Method Six.

[0154] Specific Implementation Method Eight: This implementation method further defines Specific Implementation Method Seven. In S207, the elements of the neighborhood category vector LYFCVec are as follows:

[0155] LYFCVec[1]=(MPdx-LYFCNbMeanDx) / (ImgWidth×0.1+1.0);

[0156] Among them, LYFCVec[1] is the first element of the neighborhood category vector LYFCVec;

[0157] LYFCVec[2]=(MPdy-LYFCNbMeanDy) / (ImgHeight×0.1+1.0);

[0158] Among them, LYFCVec[2] is the second element of the neighborhood category vector LYFCVec;

[0159] LYFCVec[3]=(LYFCCurMag-LYFCNbMeanMag) / (LYFCNbMeanMag+0.5);

[0160] Among them, LYFCVec[3] is the third element of the neighborhood category vector LYFCVec;

[0161] LYFCVec[4]=(LYFCCurAngle-LYFCNbMeanAngle) / π;

[0162] Among them, LYFCVec[4] is the 4th element of the neighborhood category vector LYFCVec;

[0163] If LYFCVec[4] is greater than 1.0, then let LYFCVec[4] = LYFCVec[4] - 2.0;

[0164] If LYFCVec[4] is less than -1.0, then let LYFCVec[4] = LYFCVec[4] + 2.0;

[0165] If -1.0≤LYFCVec[4]≤1.0, then LYFCVec[4]=LYFCVec[4];

[0166] LYFCVec[5]=LYFCCurMag / (LYFCMeanMag+0.5)-1.0;

[0167] Among them, LYFCVec[5] is the 5th element of the neighborhood category vector LYFCVec;

[0168] LYFCVec[6]=(LYFCCurAngle-arctan2(LYFCMeanDy,LYFCMeanDx)) / π;

[0169] Among them, LYFCVec[6] is the 6th element of the neighborhood category vector LYFCVec;

[0170] If LYFCVec[6] is greater than 1.0, then let LYFCVec[6] = LYFCVec[6] - 2.0;

[0171] If LYFCVec[6] is less than -1.0, then let LYFCVec[6] = LYFCVec[6] + 2.0;

[0172] If -1.0≤LYFCVec[6]≤1.0, then LYFCVec[6]=LYFCVec[6];

[0173]

[0174] Among them, LYFCVec[7] is the 7th element of the neighborhood category vector LYFCVec; the displacement vector of the current matching point is simultaneously projected and aligned to the horizontal direction of the UAV flight (SKBQVec[5], SKBQVec[6]) and the azimuth rotation direction (SKBQVec[1], SKBQVec[2]), which has a stronger ability to distinguish mismatched points in the edge region caused by rotation compared with using only the horizontal direction;

[0175] LYFCVec[8]=tanh(LYFCNbAngVar×5.0)×(0.5+0.5×SKBQVec[8]);

[0176] Wherein, tanh is the hyperbolic tangent function, and LYFCVec[8] is the 8th element of the neighborhood category vector LYFCVec; the formula adjusts the neighborhood displacement direction variance with the confidence weight of the time interval (between 0.5 and 1.0). The shorter the shooting time interval, the closer SKBQVec[8] is to 1, and the more reliable this dimension is. When the time interval is long, the weight is automatically reduced to avoid the inaccuracy of the direction consistency criterion under long time interval.

[0177] LYFCVec[9]=MPxA / (ImgWidth+1.0)-0.5;

[0178] Among them, LYFCVec[9] is the 9th element of the neighborhood category vector LYFCVec;

[0179] LYFCVec

[10] =MPyA / (ImgHeight+1.0)-0.5;

[0180] Among them, LYFCVec

[10] is the 10th element of the neighborhood category vector LYFCVec;

[0181] LYFCVec

[11] =(LYFCCurMag / (LYFCMeanMag×(1.0+ABS(SKBQVec[7]))+0.5))-1.0

[0182] Among them, LYFCVec

[11] is the 11th element of the neighborhood category vector LYFCVec;

[0183] LYFCVec

[11] characterizes the amplitude deviation after height correction. If SKBQVec[7]>0, the displacement amplitude of the correct matching point should also be proportionally amplified. The denominator scales the expected amplitude according to the height change. If the actual amplitude of the current point is far from the expected amplitude, the value of LYFCVec

[11] deviates from 0, and it tends to be mismatched.

[0184]

[0185] Among them, LYFCVec

[12] is the 12th element of the neighborhood category vector LYFCVec.

[0186] The other steps and parameters are the same as in Specific Implementation Method Seven.

[0187] Specific Implementation Method Nine: This implementation method is a further limitation of Specific Implementation Method Eight. The clean point set homography transformation registration module JHPYModel is used to filter all correct matching point pairs with MPIsError equal to 0 from the MatchList list, estimate the homography transformation matrix using the direct linear transformation least squares method, and then perform perspective transformation resampling on ImgB_Data to complete image registration.

[0188] The input to the clean point set homography transformation registration module JHPYModel is the list MatchList after updating the MPIsError field and the pixel data matrix ImgB_Data, where ImgB_Data is the pixel data matrix corresponding to image B;

[0189] The output of the clean point set homography transformation registration module JHPYModel is the clean registration output JHPYOutput, which specifically includes the registered image JHPYAlignImg and the homography transformation matrix JHPYTransMat.

[0190] The other steps and parameters are the same as in Specific Implementation Method 8.

[0191] Specific Implementation Method Ten: This implementation method further defines Specific Implementation Method Nine. The working process of the purification point set individual transformation registration module JHPYModel is as follows:

[0192] S301. Create a list of valid matching pairs after purification: JHPYValidList = empty list;

[0193] After initial registration, the image JHPYAlignImg is a zero-pixel matrix of the same size as ImgA_Data, and the homography transformation matrix JHPYTransMat is a 3x3 identity matrix, where ImgA_Data is the pixel data matrix corresponding to image A.

[0194] S302. Establish a valid point-to-point filtering counter JHPYCounter1=1;

[0195] S303. Retrieve the JHPYCounter1 record from the MatchList list;

[0196] If the value of the MPIsError field of the JHPYCounter1 record is equal to 0, then the matching point pair corresponding to the JHPYCounter1 record is added to the cleaned valid matching point pair list JHPYValidList.

[0197] If the value of the MPIsError field in the JHPYCounter1 record is equal to 1, then execute S304;

[0198] S304. Let JHPYCounter1 = JHPYCounter1 + 1;

[0199] If JHPYCounter1 is less than or equal to MPTotal, proceed to S303; otherwise, proceed to S305.

[0200] S305. Get the number of valid matching pairs: JHPYValidNum = the number of elements in the list JHPYValidList;

[0201] If the number of valid matching pairs JHPYValidNum is less than MinValidMatch, then increase the number of valid matching pairs manually (and add the increased valid matching pairs to the list JHPYValidList) so that the number of valid matching pairs JHPYValidNum = MinValidMatch, and then execute S306.

[0202] MinValidMatch is the minimum number of valid matching pairs, and its default value is 8.

[0203] If the number of valid matching pairs JHPYValidNum is greater than or equal to MinValidMatch, then execute S306 directly.

[0204] S306. Create a linear transformation coefficient matrix JHPYCoefMat with 2 × JHPYValidNum rows and 9 columns, and initialize all elements to 0.

[0205] Establish a linear transformation coefficient matrix filling counter JHPYCounter2=1;

[0206] S307. Retrieve the JHPYCounter2 record from the list JHPYValidList, and obtain the coordinates JHPYxA (i.e., MPxA) and JHPYyA (i.e., MPyA) of the matching point pair corresponding to the JHPYCounter2 record in image A, and the coordinates JHPYxB (i.e., MPxB) and JHPYyB (i.e., MPyB) of the matching point pair corresponding to the JHPYCounter2 record in image B.

[0207] S308. Fill in the 2×JHPYCounter2-1 row of the linear transformation coefficient matrix JHPYCoefMat. In the 2×JHPYCounter2-1 row, the values ​​of each column are -JHPYxB, -JHPYyB, -1, 0, 0, 0, JHPYxA×JHPYxB, JHPYxA×JHPYyB, and JHPYxA, respectively.

[0208] Fill in the 2nd × JHPYCounter2 row of the linear transformation coefficient matrix JHPYCoefMat, with the values ​​of each column being 0, 0, 0, -JHPYxB, -JHPYyB, -1, JHPYyA×JHPYxB, JHPYyA×JHPYyB, and JHPYyA, respectively.

[0209] S309. Let JHPYCounter2 = JHPYCounter2 + 1;

[0210] If JHPYCounter2 is less than or equal to JHPYValidNum, proceed to S307; otherwise, proceed to S310.

[0211] Where JHPYValidNum represents the total number of matching pairs in the list JHPYValidList;

[0212] S310. Perform singular value decomposition on the linear transformation coefficient matrix JHPYCoefMat to obtain the 9-dimensional right singular vector JHPYHVec corresponding to the minimum singular value; rearrange JHPYHVec into a 3x3 matrix and assign the resulting matrix to JHPYTransMat.

[0213] Specifically, the first three elements of the right singular vector JHPYHVec are taken as the first row of the homography transformation matrix, the middle three elements are taken as the second row of the homography transformation matrix, and the last three elements are taken as the third row of the homography transformation matrix.

[0214] S311, Establish perspective transformation row counter JHPYRowC=1, establish perspective transformation column counter JHPYColC=1;

[0215] S312. The homogeneous coordinate vector of the pixel in the JHPYRowC row and JHPYColC column of image A is JHPYPtA, and the three components in the homogeneous coordinate vector JHPYPtA are JHPYColC, JHPYRowC and 1, respectively.

[0216] Calculate the homogeneous coordinate vector JHPYPtB of the pixel in row JHPYRowC and column JHPYColC in image B:

[0217] JHPYPtB = the inverse matrix of JHPYTransMat multiplied by the homogeneous coordinate vector JHPYPtA

[0218] Calculate the first component of the actual column coordinate JHPYBx = JHPYPtB divided by the third component of JHPYPtB, and calculate the second component of the actual row coordinate JHPYBy = JHPYPtB divided by the third component of JHPYPtB.

[0219] S313. If 1≤JHPYBx≤ImgWidth and 1≤JHPYBy≤ImgHeight are satisfied, then extract the pixel value of the JHPYBy row and JHPYBx column in the pixel data matrix ImgB_Data (using bilinear interpolation to obtain the pixel at the decimal coordinates), and assign the extracted pixel value to the pixel of the JHPYRowC row and JHPYColC column in the registered image JHPYAlignImg.

[0220] If 1≤JHPYBx≤ImgWidth and 1≤JHPYBy≤ImgHeight are not satisfied, then keep the pixel value of the JHPYRowC column in the registered image JHPYAlignImg as 0.

[0221] S314. Let JHPYColC = JHPYColC + 1;

[0222] If JHPYColC is less than or equal to ImgWidth, then proceed to S312;

[0223] Otherwise, set JHPYColC=1, and let JHPYRowC=JHPYRowC+1. If JHPYRowC is less than or equal to ImgHeight, then go to S312; otherwise, go to S315.

[0224] S315, Output the registered image JHPYAlignImg, S3 ends.

[0225] The other steps and parameters are the same as in Specific Implementation Method Nine.

[0226] The above examples of the present invention are merely illustrative of the computational model and process of the present invention, and are not intended to limit the implementation of the present invention. Those skilled in the art will recognize that other variations or modifications can be made based on the above description. It is impossible to exhaustively list all possible implementations here. Any obvious variations or modifications derived from the technical solutions of the present invention are still within the scope of protection of the present invention.

Claims

1. A UAV image matching method based on spatiotemporal changes and neighborhood category outliers, characterized in that, The method specifically includes the following steps: Step 1: Acquire images A and B to be registered, and record the flight parameters of the UAV when capturing images A and B respectively. The flight parameters include the azimuth angle, pitch angle, flight altitude, X coordinate of the UAV center, Y coordinate of the UAV center, and capture time when capturing the images. Perform initial matching on image A and image B to obtain an initial matching point pair list MatchList. Each record in the MatchList includes the following fields: The MPID field is the number of the matching point pair; The MPxA field represents the pixel coordinates of the matched point pair in the width direction within image A. The MPyA field represents the pixel coordinates of the matching point pair in image A in the height direction. The MPxB field represents the pixel coordinates of the matched point pair in the width direction within image B. The MPyB field represents the pixel coordinates of the matching point pair in the height direction within image B. The MPdx field represents the value of MPxB - MPxA; The MPdy field represents the value of MPyB-MPyA; And record the pixel data matrices corresponding to image A and image B respectively; Step 2: Based on the initial matching point pair list MatchList, the pixel data matrix corresponding to image A and image B, and using the spatiotemporal change description vector construction module SKBQModel, the neighborhood category outlier matching point detection module LYFCModel, and the clean point set homography transformation registration module JHPYModel, the registered image is obtained. The input to the spatiotemporal change description vector construction module SKBQModel is the flight parameters when the UAV captures images A and B, and the output is the spatiotemporal change description vector SKBQVec. The neighborhood category outlier matching point detection module LYFCModel takes the initial matching point pair list MatchList and the spatiotemporal change description vector SKBQVec as input and outputs the error point pair list LYFCErrorList. The input to the clean point set homography transformation registration module JHPYModel is the list MatchList after updating the MPIsError field and the pixel data matrix ImgB_Data, where ImgB_Data is the pixel data matrix corresponding to image B; The output of the clean point set homography transformation registration module JHPYModel is the clean registration output JHPYOutput, which specifically includes the registered image JHPYAlignImg and the homography transformation matrix JHPYTransMat.

2. The UAV image matching method based on spatiotemporal change and neighborhood category outlier as described in claim 1, characterized in that, The specific working process of the spatiotemporal transition description vector construction module SKBQModel is as follows: S101. Initialize the spatiotemporal transition description vector SKBQVec = an 8-dimensional all-zero vector; S102, Calculate the azimuth difference in radians SKBQTemp1: SKBQTemp1=(ImgB_Azi-ImgA_Azi)×π / 180 Where π is the mathematical constant pi, ImgA_Azi is the azimuth angle when the drone captures image A, and ImgB_Azi is the azimuth angle when the drone captures image B. Calculate the pitch angle difference in radians using SKBQTemp2: SKBQTemp2=(ImgB_Ele-ImgA_Ele)×π / 180 Where ImgA_Ele is the pitch angle when the drone captures image A, and ImgB_Ele is the pitch angle when the drone captures image B; S103, let SKBQVec[1]=sin(SKBQTemp1); Wherein, SKBQVec[1] is the first element in the spatiotemporal change description vector SKBQVec, and sin is the sine function; Let SKBQVec[2] = cos(SKBQTemp1); Wherein, cos is the cosine function, and SKBQVec[2] is the second element in the spatiotemporal change description vector SKBQVec; Let SKBQVec[3]=sin(SKBQTemp2), SKBQVec[4]=cos(SKBQTemp2); Among them, SKBQVec[3] is the third element in the spatiotemporal change description vector SKBQVec, and SKBQVec[4] is the fourth element in the spatiotemporal change description vector SKBQVec. S104. Calculate the normalized horizontal displacement SKBQTemp3 and the normalized vertical displacement SKBQTemp4: Where ImgA_CX is the X coordinate of the drone's center when image A is captured, ImgB_CX is the X coordinate of the drone's center when image B is captured, ABS is the absolute value, ImgA_CY is the Y coordinate of the drone's center when image A is captured, and ImgB_CY is the Y coordinate of the drone's center when image B is captured. Let SKBQVec'[5]=SKBQTemp3, SKBQVec[6]=SKBQTemp4; Among them, SKBQVec'[5] is the 5th element in the spatiotemporal change description vector SKBQVec, and SKBQVec[6] is the 6th element in the spatiotemporal change description vector SKBQVec; S105, Calculate the normalized value of height change SKBQTemp5: SKBQTemp5=(ImgB_Alt-ImgA_Alt) / (ABS(ImgB_Alt-ImgA_Alt)+5.0) Where ImgA_Alt is the altitude when the drone captured image A, and ImgB_Alt is the altitude when the drone captured image B; Let SKBQVec[7] = SKBQTemp5; Among them, SKBQVec[7] is the 7th element in the spatiotemporal change description vector SKBQVec; S106, Calculate the time interval attenuation SKBQTemp6: SKBQTemp6=ABS(ImgB_Time-ImgA_Time) Where ImgA_Time is the time when image A was captured, and ImgB_Time is the time when image B was captured; Let SKBQVec[8]=exp(-SKBQTemp6 / TimeTau); Where exp represents an exponential function with the natural constant e as the base, TimeTau is the time decay constant, and SKBQVec[8] is the 8th element in the spatiotemporal change description vector SKBQVec; S107, Calculate the azimuth-elevation coupling cross term SKBQTemp7: SKBQTemp7=sin(SKBQTemp1)×cos(SKBQTemp2)×0.5 The 5th element in the spacetime transition description vector SKBQVec is updated as follows: SKBQVec[5]=SKBQVec'[5]+SKBQTemp7×0.3 Wherein, SKBQVec[5] represents the 5th element in the final spatiotemporal transition description vector SKBQVec; S108. Output the spatiotemporal transition description vector SKBQVec as the result of the SKBQModel module.

3. The UAV image matching method based on spatiotemporal change and neighborhood category outlier as described in claim 2, characterized in that, The working process of the neighborhood category outlier matching point detection module LYFCModel is as follows: S201. Initialize the error point pair list LYFCErrorList = empty list; S202. Calculate the global mean of the displacement vector in the width direction and the global mean of the displacement vector in the height direction, and then calculate the global displacement mean amplitude LYFCMeanMag based on the global mean of the displacement vector in the width direction and the global mean of the displacement vector in the height direction. S203. Establish a point-by-point category vector calculation counter LYFCCounter1=1; S204. Retrieve the LYFCCounter1 record from the MatchList and obtain the MPxA, MPyA, MPdx, and MPdy fields of the LYFCCounter1 record. Calculate the displacement magnitude LYFCCurMag of the matching point pair corresponding to the LYFCCounter1 record: Calculate the displacement direction LYFCCurAngle of the matching point pair corresponding to the LYFCCounter1 record: Where arctan2 is the arctangent function in the four quadrants; S205. Calculate the Euclidean distance between each other matching point in image A and the matching point corresponding to the LYFCCounter1 record in image A, and store the NbK matching points with the smallest distances into the neighborhood list LYFCNbList. Where NbK is the number of nearest neighbors in the neighborhood; S206. Calculate the mean displacement of the neighborhood in the width direction (LYFCNbMeanDx) and the mean displacement in the height direction (LYFCNbMeanDy) of the corresponding matching point of the LYFCCounter1 record in image A, and then calculate the magnitude of the mean displacement of the neighborhood (LYFCNbMeanMag) of the corresponding matching point of the LYFCCounter1 record in image A. And calculate the mean of the neighborhood displacement direction LYFCNbMeanAngle and the variance of the neighborhood displacement direction LYFCNbAngVar of the corresponding matching point in image A for the LYFCCounter1 record; S207. Construct the 12-dimensional neighborhood category vector LYFCVec for the matching point pair corresponding to the LYFCCounter1 record in the list MatchList; S208. Assign the neighborhood category vector LYFCVec to the MPNbVec field of the LYFCCounter1 record in the list MatchList. Let LYFCCounter1 = LYFCCounter1 + 1; If LYFCCounter1 is less than or equal to MPTotal, proceed to S204; otherwise, proceed to S209. S209. Set the number of clusters to ClusterNum, perform K-means clustering on the MPNbVec field of all records in the list MatchList, obtain the cluster label MPCluster for each matching pair based on the clustering results, and obtain the matrix LYFCCenters composed of the cluster centers of each cluster. S210. Count the number of matching point pairs contained in each cluster, select the cluster with the most matching point pairs as the dominant cluster DominantCluster, and denote the cluster center vector corresponding to the dominant cluster as DominantCenter. S211. Establish an outlier detection counter LYFCCounter3=1; S212. Retrieve the LYFCCounter3rd record from the MatchList list; If the clustering label MPCluster of the matching point pair corresponding to the LYFCCounter3 record is equal to DominantCluster, then proceed to S213; If the clustering label MPCluster of the matching point pair corresponding to the LYFCCounter3 record is not equal to DominantCluster, then calculate the Euclidean distance LYFCDistToDom between the MPNbVec field of the matching point pair corresponding to the LYFCCounter3 record and the DominantCenter. The distance from the MPNbVec field to the DominantCenter for each matching point pair within the dominant cluster is calculated, and the standard deviation LYFCDomStd of the calculated distance is then calculated. If LYFCDistToDom is greater than LYFCDomStd×OutlierMult, then set the error flag MPIsError=1 for the matching point pair corresponding to the LYFCCounter3 record; add the matching point pair corresponding to the LYFCCounter3 record to the error point pair list LYFCErrorList; If LYFCDistToDom is less than or equal to LYFCDomStd×OutlierMult, then set the error flag MPIsError=0 for the matching point pair corresponding to the LYFCCounter3 record; Where OutlierMult is the outlier distance factor; S213. Set the error flag MPIsError=0 for the matching point pair corresponding to the LYFCCounter3 record, and let LYFCCounter3=LYFCCounter3+1; If LYFCCounter3 is less than or equal to MPTotal, go to S212; otherwise, go to S214. S214, S2 ends.

4. The UAV image matching method based on spatiotemporal change and neighborhood category outlier as described in claim 3, characterized in that, The specific process of S202 is as follows: Step 1: Establish a global mean temporary storage value for the width direction displacement vector LYFCMeanDx'=0 and a global mean temporary storage value for the height direction displacement vector LYFCMeanDy'=0; And establish a summation counter LYFCCounter0=1; Step 2: Set LYFCMeanDx' = LYFCMeanDx' + the value of the MPdx field of the LYFCCounter0th record in the list MatchList; Let LYFCMeanDy' = LYFCMeanDy' + the value of the MPdy field of the LYFCCounter0th record in the list MatchList; Step 3: Set LYFCCounter0 = LYFCCounter0 + 1; If LYFCCounter0 is less than or equal to the total number of records in the MatchList (MPTotal), proceed to step 2; otherwise, proceed to step 4. Step 4: Final global mean of the width direction displacement vector LYFCMeanDx: LYFCMeanDx=LYFCMeanDx' / MPTotal The final global mean of the height displacement vector is LYFCMeanDy. LYFCMeanDy=LYFCMeanDy' / MPTotal Step 5: Calculate the global mean displacement magnitude LYFCMeanMag based on LYFCMeanDx and LYFCMeanDy: 。 5. The UAV image matching method based on spatiotemporal change and neighborhood category outlier as described in claim 4, characterized in that, The specific process of S206 is as follows: Step 1) Initialize the temporary storage of the average displacement of the neighborhood in the width direction LYFCNbMeanDx'=0, the temporary storage of the average displacement of the neighborhood in the height direction LYFCNbMeanDy'=0, and the neighborhood displacement direction list LYFCNbAngles=empty list; Set the neighborhood traversal counter LYFCCounter2=1; Step 2) Extract the MPdx and MPdy fields corresponding to the LYFCCounter2 neighbor point in the list LYFCNbList in the list MatchList; Let LYFCNbMeanDx'=LYFCNbMeanDx'+ ; LYFCNbMeanDy’=LYFCNbMeanDy’+ ; Calculate the displacement direction of the LYFCCounter2th neighboring point in the list LYFCNbList. ,Will Add to list LYFCNbAngles; in, The value of the MPdx field for the LYFCCounter2 neighboring point. The value of the MPdy field for the LYFCCounter2 neighboring point. The displacement direction of the LYFCCounter2th neighboring point; Step 3) Set LYFCCounter2 = LYFCCounter2 + 1; If LYFCCounter2 is less than or equal to NbK, proceed to step 2; otherwise, proceed to step 4. Step 4) Calculate the mean displacement of the neighborhood in the width direction, LYFCNbMeanDx: LYFCNbMeanDx=LYFCNbMeanDx' / NbK Calculate the mean displacement of the neighborhood in the height direction, LYFCNbMeanDy: LYFCNbMeanDy=LYFCNbMeanDy' / NbK Calculate the mean displacement magnitude of the neighborhood LYFCNbMeanMag: Calculate the mean displacement direction of the neighborhood LYFCNbMeanAngle: Step 5) Calculate the neighborhood displacement direction variance LYFCNbAngVar: in, This refers to the first element in the list LYFCNbAngles. Values.

6. The UAV image matching method based on spatiotemporal change and neighborhood category outlier as described in claim 5, characterized in that, In S207, the elements of the neighborhood category vector LYFCVec are as follows: LYFCVec[1]=(MPdx-LYFCNbMeanDx) / (ImgWidth×0.1+1.0); Among them, LYFCVec[1] is the first element of the neighborhood category vector LYFCVec; LYFCVec[2]=(MPdy-LYFCNbMeanDy) / (ImgHeight×0.1+1.0); Among them, LYFCVec[2] is the second element of the neighborhood category vector LYFCVec; LYFCVec[3]=(LYFCCurMag-LYFCNbMeanMag) / (LYFCNbMeanMag+0.5); Among them, LYFCVec[3] is the third element of the neighborhood category vector LYFCVec; LYFCVec[4]=(LYFCCurAngle-LYFCNbMeanAngle) / π; Among them, LYFCVec[4] is the 4th element of the neighborhood category vector LYFCVec; If LYFCVec[4] is greater than 1.0, then let LYFCVec[4] = LYFCVec[4] - 2.0; If LYFCVec[4] is less than -1.0, then let LYFCVec[4] = LYFCVec[4] + 2.0; If -1.0≤LYFCVec[4]≤1.0, then LYFCVec[4]=LYFCVec[4]; LYFCVec[5]=LYFCCurMag / (LYFCMeanMag+0.5)-1.0; Among them, LYFCVec[5] is the 5th element of the neighborhood category vector LYFCVec; LYFCVec[6]=(LYFCCurAngle-arctan2(LYFCMeanDy,LYFCMeanDx)) / π; Among them, LYFCVec[6] is the 6th element of the neighborhood category vector LYFCVec; If LYFCVec[6] is greater than 1.0, then let LYFCVec[6] = LYFCVec[6] - 2.0; If LYFCVec[6] is less than -1.0, then let LYFCVec[6] = LYFCVec[6] + 2.0; If -1.0≤LYFCVec[6]≤1.0, then LYFCVec[6]=LYFCVec[6]; Among them, LYFCVec[7] is the 7th element of the neighborhood category vector LYFCVec; LYFCVec[8]=tanh(LYFCNbAngVar×5.0)×(0.5+0.5×SKBQVec[8]); Where tanh is the hyperbolic tangent function, and LYFCVec[8] is the 8th element of the neighborhood category vector LYFCVec; LYFCVec[9]=MPxA / (ImgWidth+1.0)-0.5; Among them, LYFCVec[9] is the 9th element of the neighborhood category vector LYFCVec; LYFCVec[10]=MPyA / (ImgHeight+1.0)-0.5; Among them, LYFCVec[10] is the 10th element of the neighborhood category vector LYFCVec; LYFCVec[11]=(LYFCCurMag / (LYFCMeanMag×(1.0+ABS(SKBQVec[7]))+0.5))-1.0 Among them, LYFCVec[11] is the 11th element of the neighborhood category vector LYFCVec; Among them, LYFCVec[12] is the 12th element of the neighborhood category vector LYFCVec.

7. The UAV image matching method based on spatiotemporal change and neighborhood category outlier as described in claim 6, characterized in that, The working process of the purification point set single-application transformation registration module JHPYModel is as follows: S301. Create a list of valid matching pairs after purification: JHPYValidList = empty list; After initial registration, the image JHPYAlignImg is a zero-pixel matrix of the same size as ImgA_Data, and the homography transformation matrix JHPYTransMat is a 3x3 identity matrix, where ImgA_Data is the pixel data matrix corresponding to image A. S302. Establish a valid point-to-point filtering counter JHPYCounter1=1; S303. Retrieve the JHPYCounter1 record from the MatchList list; If the value of the MPIsError field of the JHPYCounter1 record is equal to 0, then the matching point pair corresponding to the JHPYCounter1 record is added to the cleaned valid matching point pair list JHPYValidList. If the value of the MPIsError field in the JHPYCounter1 record is equal to 1, then execute S304; S304. Let JHPYCounter1 = JHPYCounter1 + 1; If JHPYCounter1 is less than or equal to MPTotal, proceed to S303; otherwise, proceed to S305. S305. Get the number of valid matching pairs: JHPYValidNum = the number of elements in the list JHPYValidList; If the number of valid matching pairs JHPYValidNum is less than MinValidMatch, then increase the number of valid matching pairs manually to make the number of valid matching pairs JHPYValidNum = MinValidMatch, and then execute S306. Wherein, MinValidMatch is the minimum number of valid matching pairs; If the number of valid matching pairs JHPYValidNum is greater than or equal to MinValidMatch, then execute S306 directly. S306. Create a linear transformation coefficient matrix JHPYCoefMat with 2 × JHPYValidNum rows and 9 columns, and initialize all elements to 0. Establish a linear transformation coefficient matrix filling counter JHPYCounter2=1; S307. Retrieve the JHPYCounter2 record from the list JHPYValidList, obtain the coordinates JHPYxA and JHPYyA of the matching point pair corresponding to the JHPYCounter2 record in image A, and the coordinates JHPYxB and JHPYyB of the matching point pair corresponding to the JHPYCounter2 record in image B. S308. Fill in the 2×JHPYCounter2-1 row of the linear transformation coefficient matrix JHPYCoefMat. In the 2×JHPYCounter2-1 row, the values ​​of each column are -JHPYxB, -JHPYyB, -1, 0, 0, 0, JHPYxA×JHPYxB, JHPYxA×JHPYyB, and JHPYxA, respectively. Fill in the 2nd × JHPYCounter2 row of the linear transformation coefficient matrix JHPYCoefMat, with the values ​​of each column being 0, 0, 0, -JHPYxB, -JHPYyB, -1, JHPYyA×JHPYxB, JHPYyA×JHPYyB, and JHPYyA, respectively. S309. Let JHPYCounter2 = JHPYCounter2 + 1; If JHPYCounter2 is less than or equal to JHPYValidNum, proceed to S307; otherwise, proceed to S310. Where JHPYValidNum represents the total number of matching pairs in the list JHPYValidList; S310. Perform singular value decomposition on the linear transformation coefficient matrix JHPYCoefMat to obtain the 9-dimensional right singular vector JHPYHVec corresponding to the minimum singular value; rearrange JHPYHVec into a 3x3 matrix and assign the resulting matrix to JHPYTransMat. S311, Establish perspective transformation row counter JHPYRowC=1, establish perspective transformation column counter JHPYColC=1; S312. The homogeneous coordinate vector of the pixel in the JHPYRowC row and JHPYColC column of image A is JHPYPtA, and the three components in the homogeneous coordinate vector JHPYPtA are JHPYColC, JHPYRowC and 1, respectively. Calculate the homogeneous coordinate vector JHPYPtB of the pixel in row JHPYRowC and column JHPYColC in image B: JHPYPtB = the inverse matrix of JHPYTransMat multiplied by the homogeneous coordinate vector JHPYPtA Calculate the first component of the actual column coordinate JHPYBx = JHPYPtB divided by the third component of JHPYPtB, and calculate the second component of the actual row coordinate JHPYBy = JHPYPtB divided by the third component of JHPYPtB. S313. If 1≤JHPYBx≤ImgWidth and 1≤JHPYBy≤ImgHeight are satisfied, then the pixel value in the JHPYBy row and JHPYBx column of the pixel data matrix ImgB_Data is extracted and the extracted pixel value is assigned to the pixel in the JHPYRowC row and JHPYColC column of the registered image JHPYAlignImg. If 1≤JHPYBx≤ImgWidth and 1≤JHPYBy≤ImgHeight are not satisfied, then keep the pixel value of the JHPYRowC column in the registered image JHPYAlignImg as 0. S314. Let JHPYColC = JHPYColC + 1; If JHPYColC is less than or equal to ImgWidth, then proceed to S312; Otherwise, set JHPYColC=1, and let JHPYRowC=JHPYRowC+1. If JHPYRowC is less than or equal to ImgHeight, then go to S312; otherwise, go to S315. S315, Output the registered image JHPYAlignImg, S3 ends.

Images