A method for registering large inclination aerial images and orthophotos

Abstract

The present invention relates to a method for registering large tilt aerial images and orthophotos. The technical solution is as follows: First, use the ASIFT algorithm to simulate the affine image perspective of the large tilt aerial image S to obtain the simulated affine image S'; then construct a non-linear scale space for the simulated affine image S' and the orthophoto R; in the non-linear scale space, use phase consistency and the Sobel operator to construct a feature point metric matrix for feature point extraction, and construct a GLOH-like feature description operator by calculating the amplitude and direction of the feature points; then perform feature registration and error elimination on the simulated affine image S' and the orthophoto R to obtain correctly registered corresponding points. The corresponding points in the large tilt aerial image S are obtained from the correctly registered corresponding points of the affine simulated image S' through the simulation parameters, and finally, the large tilt aerial image S and the orthophoto R are corrected according to the affine transformation parameters Ψ. The present invention has the characteristics of small registration error, good anti-noise performance, high success rate and stable results.

Description

Technical Field

[0001] The present invention belongs to the technical field of registration of tilt images and orthophoto images. In particular, it relates to a method for registering large-tilt aerial images and orthophoto images. Background Art

[0002] Compared with orthophoto images, tilt images can observe ground objects from multiple angles, more truly reflect the actual situation of ground objects, and greatly compensate for the deficiencies of applications based on orthophoto images. The fusion application of tilt images and orthophoto images can quickly associate, integrate and analyze ground object information, which is of great significance to the development of remote sensing interpretation technology. As the first step in remote sensing image data fusion, image registration lays a foundation for remote sensing image earth observation applications. Image registration refers to the process of geometric alignment and registration of two or more images of the same scene obtained by different sensors, different perspectives and different times. The accuracy of image registration plays a crucial role in subsequent applications. How to achieve high-efficiency and high-precision registration of large-tilt aerial images is an issue of concern to those skilled in the art.

[0003] Image registration can generally be divided into region-based methods and feature-based methods. Region-based methods rely on georeferencing technology to roughly register images and eliminate obvious translation and rotation differences between image pairs. However, this method is not suitable for large-tilt aerial images because there are very serious geometric deformations of ground objects in large-tilt aerial images, and it is difficult to obtain their geographic coordinate information. Feature-based methods achieve image registration by evaluating the structural feature information of images. This method is more controllable and robust. At present, common registration methods are generally applicable to the registration of images with small tilt angles, but are not applicable to the registration of large-tilt aerial images. There are relatively few patented technologies and papers on the registration of large-tilt aerial images.

[0004] The literature 1 (Lowe D G. Distinctive image features from scale-invariant keypoints[J]. International Journal of Computer Vision, 2004, 60: 91-110.) is the Scale-Invariant Feature Transform (SIFT). This method can handle image matching well under scale and rotation conditions. However, for feature point matching of oblique images, only a small number of features can be extracted, and the registration error is large. The literature 2 (Zhao X, Zhu Q, Xiao X, et al. Automatic matching method for aviation oblique images based on homography transformation[J]. Journal of Computer Applications, 2015, 35(6): 1720.) proposed an automatic registration method for large-angle aviation oblique images, H-SIFT. The homography transformation matrix between two images A and B is calculated using the rough exterior orientation elements of the images. The corrected image A' of the original image A is obtained using the homography transformation matrix. Then, SIFT matching is performed on the corrected image A' and the image B, and the feature points matched on the corrected image A' are inversely calculated to the original image A. This method highly depends on the distribution and accuracy of the initial matching points, and the matching result is unstable. The literature 3 (Morel J M, Yu G. ASIFT: A new framework for fully affine invariant image comparison[J]. SIAM Journal on Imaging Sciences, 2009, 2(2): 438-469.) proposed an affine-invariant image feature matching algorithm (ASIFT), which solved the matching problem of SIFT for oblique images. ASIFT has complete affine invariance and achieved good results in the registration of large-angle images. However, the registration accuracy of this method for large-angle aviation images is low. The literature 4 (Zhao Chaohe, Yang Huachao, Zhang Lei, et al. An automatic image registration algorithm for large angles at sub-pixel level[J]. Bulletin of Surveying and Mapping, 2014(8): 30-35.) utilized the complete radiometric invariance of ASIFT and proposed a sub-pixel automatic registration algorithm for large-angle images based on the integrated complementary invariant features of ASIFT and Harris. However, there are serious problems in large-angle aviation images, such as severe geometric deformation of ground objects, resolution differences, image rotation, and mutual occlusion of ground objects. The methods based on SIFT and ASIFT use Gaussian kernels to construct the scale space and rely on the construction of feature description operators using image gray information and gradient features, and cannot handle well the problem of complex ground object information in aviation images and extract rich structural features.In addition, the gray-scale information and gradient features are sensitive to image noise, and the constructed description operator is not robust, resulting in a low registration success rate and large errors. Summary of the Invention

[0005] The present invention aims to overcome the defects of the prior art, and the purpose is to provide a method for registering large-angle aerial images and orthophotos with small registration errors, good anti-noise performance, high success rate and stable results.

[0006] In order to achieve the above purpose, the technical solution adopted by the present invention is:

[0007] Step 1: Use the ASIFT algorithm to construct a simulated affine image S' of the large-angle aerial image S

[0008] Step 1.1: Parameter sampling

[0009] According to the ASIFT algorithm, in a hemispherical space, the viewing angle of the large-angle aerial image S is simulated through the direction parameters of the camera in the longitude direction angle φ and the latitude direction angle θ; the longitude direction angle φ and the latitude direction angle θ are discretized respectively. Among them:

[0010] The latitude direction angle θ adopts a geometric series discretization method:

[0011] T n = a n In formula (1):

[0012] T n represents the inclination of the camera at the nth sampling level in the latitude direction;

[0013] n represents the sampling level of the camera in the latitude direction, n = 1,..., N;

[0014] a represents the sampling interval of the camera in the latitude direction, a > 1.

[0015] The longitude direction angle φ adopts an arithmetic series discretization method:

[0016]

[0017] In formula (2):

[0018] represents the rotation angle at the kth sampling level in the longitude direction under the inclination T of the camera in the latitude direction; n

[0019] b / T n represents the sampling interval in the longitude direction under the inclination T of the camera in the latitude direction, and b is an angle constant; n

[0020] ​​k represents the tilt angle T of the camera in the latitude direction n Under this condition, the sampling level in the longitude direction, k = 1,..., K

[0021]

[0022] Step 1.2: Construct the simulated affine image S' of the large tilt angle aerial image S

[0023] Let the tilt angle of the camera in the latitude direction be T n , and the rotation angle in the longitude direction be The construction process of the simulated affine image of the large tilt angle aerial image S is as follows: First, perform a bilinear interpolation rotation operation on the large tilt angle aerial image S to obtain a rotated image; then, perform convolution on the rotated image with an anti-aliasing filter to obtain the rotated image after filtering; then, perform a tilt operation on the rotated image after filtering to obtain the simulated affine image of the large tilt angle aerial image S when the tilt angle of the camera in the latitude direction is T n and the rotation angle in the longitude direction is the simulated affine image of the large tilt angle aerial image S

[0024]

[0025] In formula (3):

[0026] represents the simulated affine image of the large tilt angle aerial image S when the tilt angle of the camera in the latitude direction is T n and the rotation angle in the longitude direction is ;

[0027] Tilt represents performing a tilt operation on the rotated image after filtering with an anti-aliasing filter;

[0028] G δ represents the anti-aliasing filter, where δ represents the kernel standard deviation of the anti-aliasing filter, δ = 0.8T n ;

[0029] * represents the convolution operation;

[0030] Rot represents performing a bilinear interpolation rotation operation on the large tilt angle aerial image S;

[0031] S represents the large tilt angle aerial image;

[0032] represents the tilt angle T of the camera in the latitude direction n Under this condition, the rotation angle at the k-th sampling level in the longitude direction;

[0033] T nIndicates the inclination of the camera at the nth sampling level in the latitude direction.

[0034] And so on, the remaining simulated affine images of the large-inclination aerial image S are obtained. Then, the simulated affine image

[0035] For simplicity of description: hereinafter, the "simulated affine image S' of the large-inclination aerial image S" will be abbreviated as "simulated affine image S'"; the "inclination of the camera in the latitude direction is T n and the rotation angle in the longitude direction is when the simulated affine image of the large-inclination aerial image S " will be abbreviated as "sub-simulated affine image ".

[0036] Step 2: Register the simulated affine image S' and the orthoimage R

[0037] Step 2.1: Construct the non-linear scale space

[0038] Discretize the non-linear scale space of the sub-simulated affine image into 1 base space and L sub-spaces. After discretization, the scales corresponding to the non-linear scale space of the sub-simulated affine image are as follows:

[0039]

[0040] In formula (4):

[0041] i represents the i-th layer of the non-linear scale space of the sub-simulated affine image , i = 0, 1, 2,...., L;

[0042] σ i represents the scale of the i-th layer non-linear scale space of the sub-simulated affine image ;

[0043] σ0 represents the initial scale of the non-linear scale space of the sub-simulated affine image , σ0 is a constant;

[0044] L represents the number of sub-spaces of the sub-simulated affine image .

[0045] Let the non-linear scale space of the sub-simulated affine image be: Where: represents the base space of the sub-simulated affine image , successively represent the sub-simulated affine image The first sub - space, ……, the L - th sub - space.

[0046] Sub - simulated affine image The base space of The construction process is as follows:

[0047] Use a Gaussian filter to smooth the sub - simulated affine image The Gaussian filter has a Gaussian kernel standard deviation that is the same as the initial scale σ0 of the non - linear scale space of the sub - simulated affine image to obtain the base space of the sub - simulated affine image The base space of

[0048] Sub - simulated affine image The j - th sub - space of The construction process is as follows:

[0049] First, down - sample the (j - 1) - th layer of the non - linear scale space of the sub - simulated affine image at a sampling rate of 1 / σ0 to obtain a sampling space Then, use anisotropic diffusion filtering to filter the sampling space to obtain the j - th sub - space of the sub - simulated affine image The j - th sub - space of

[0050]

[0051] In Equation (5):

[0052] represents the anisotropic diffusion filtering space obtained under the condition of the time step ;

[0053] j represents the j - th of the sub - spaces of the sub - simulated affine image , where j = 1, 2, …, L;

[0054] ↓ represents the down - sampling operation;

[0055] represents the sampling space obtained by down - sampling the (j - 1) - th layer of the non - linear scale space of the sub - simulated affine image ;

[0056] t represents the time measure;

[0057] div represents the divergence operator;

[0058] D(x, y, t) represents the diffusion coefficient, and the diffusion coefficient depends on the gradient norm

[0059] denotes the gradient operator;

[0060] denotes the j-th child space of the sub-simulation affine image ;

[0061] I d denotes the identity matrix;

[0062] σ j denotes the scale of the j-th layer of the non-linear scale space of the sub-simulation affine image ;

[0063] σ j-1 denotes the scale of the (j - 1)-th layer of the non-linear scale space of the sub-simulation affine image ;

[0064] l denotes the direction axis of anisotropic diffusion;

[0065] m denotes the number of direction axes of anisotropic diffusion;

[0066] A l denotes the diffusion coefficient matrix of the sampling space along the direction axis l of anisotropic diffusion, and A l is the discrete form of the diffusion coefficient D.

[0067] According to Equation (5), the L child spaces of the sub-simulation affine image are obtained; then the non-linear scale space of the sub-simulation affine image is: And so on, the non-linear scale spaces of the remaining sub-simulation affine images in the simulated affine image S' are obtained.

[0068] Referring to the method for obtaining the non-linear scale spaces of the remaining sub-simulation affine images in the simulated affine image S', the non-linear scale space of the orthoimage R is: R0, R1,..., R L .

[0069] Step 2.2, Feature point detection

[0070] For the sake of convenience of description, the phase consistency is denoted by "PC".

[0071] Feature point detection is performed on the sub-simulation affine image , and the feature point detection process is as follows:

[0072] First, use PC to extract the structural features of the image of the i-th layer of the non-linear scale space of the sub-simulation affine image to obtain the PC feature map of the i-th layer of the non-linear scale space of the sub-simulation affine image ; then use the Sobel operator templates in the X direction and the Y direction for the sub-simulation affine image Convolve the PC feature map of the i-th layer of the non-linear scale space to obtain a sub-simulated affine image The first-order derivatives of the PC feature map of the i-th layer of the non-linear scale space in the X and Y directions:

[0073]

[0074]

[0075] In equations (6) to (7):

[0076] i represents the i-th layer of the non-linear scale space of the sub-simulated affine image where i = 0, 1, 2,..., L;

[0077] represents the sub-simulated affine image The first-order derivative of the PC feature map of the i-th layer of the non-linear scale space in the X direction;

[0078] represents the sub-simulated affine image The first-order derivative of the PC feature map of the i-th layer of the non-linear scale space in the Y direction;

[0079] PC i represents the sub-simulated affine image The PC feature map of the i-th layer of the non-linear scale space;

[0080] * represents the convolution operation;[[ID=4o]]

[0081] Γ x represents the Sobel operator template in the X direction;

[0082] Γ y represents the Sobel operator template in the Y direction.

[0083] Then, use the first-order derivatives of the PC feature map of the i-th layer of the non-linear scale space of the obtained sub-simulated affine image in the X and Y directions to construct the feature point metric matrix M of the i-th layer of the non-linear scale space of the sub-simulated affine image : i :

[0084]

[0085] In equation (8):

[0086] i represents the i-th layer of the non-linear scale space of the sub-simulated affine image where i = 0, 1, 2,...., L;

[0087] Mi Denote the feature point metric matrix of the \(i\)-th layer of the non-linear scale space of the sub-simulated affine image ;

[0088] Denote the Gaussian kernel, and the standard deviation of the Gaussian kernel is where \(\sigma\) i Denote the scale of the \(i\)-th layer of the scale space of the sub-simulated affine image ;

[0089] * denotes the convolution operation;

[0090] Denote the first-order derivative in the \(X\) direction of the PC feature map of the \(i\)-th layer of the non-linear scale space of the sub-simulated affine image ;

[0091] Denote the first-order derivative in the \(Y\) direction of the PC feature map of the \(i\)-th layer of the non-linear scale space of the sub-simulated affine image ;

[0092] Finally, perform eigen-decomposition on the feature point metric matrix \(M\) of the \(i\)-th layer of the non-linear scale space of the sub-simulated affine image to obtain the first eigenvalue \(\lambda_1\) and the second eigenvalue \(\lambda_2\) of the feature point metric matrix \(M\) of the \(i\)-th layer of the non-linear scale space of the sub-simulated affine image i The feature point intensity map of the \(i\)-th layer of the non-linear scale space of the sub-simulated affine image is: \(R\) ; i \(=\min(\lambda_1,\lambda_2)\), and the points satisfying \(R\) \(>r\) are the feature points of the \(i\)-th layer of the non-linear scale space of the sub-simulated affine image i , where \(r\) is a threshold constant i ;

[0093] Among the feature points of each layer of the non-linear scale space of the obtained sub-simulated affine image , if there are more than 2 feature points at the same position in the non-linear scale spaces of different layers, only one of these feature points is retained; if there are feature points at different positions in the non-linear scale spaces of different layers, all of these feature points are retained, and the feature point set of the sub-simulated affine image is obtained

[0094] And so on, perform feature point detection on the remaining sub-simulated affine images in the simulated affine image \(S'\).

[0095] Combine the feature point sets of all sub-simulated affine images to obtain the feature point set of the simulated affine image \(S'\) ​

[0096] Obtain the feature point set Key of the orthophoto R by referring to the method for obtaining the feature point set of the simulated affine image S'. R ={(x1′,y1′),...,(x' m ,y' m )}

[0097] Step 2.3, Construction of feature description operator

[0098] For any feature point (x c ,y c ) in the simulated affine image S', extract the amplitude Mag c and direction Ang c of the circular region Ω centered on the feature point (x Ω and direction Ang Ω :

[0099]

[0100]

[0101] In formulas (9) to (10): <(

[0102] Ω represents a circular region centered on the feature point (x c ,y c ), and the radius of the circular region Ω is d r ;

[0103] Mag Ω represents the amplitude of the circular region Ω centered on the feature point (x c ,y c );

[0104] represents the first derivative of the PC feature map of the circular region Ω centered on the feature point (x c ,y c ) in the X direction;

[0105] represents the first derivative of the PC feature map of the circular region Ω centered on the feature point (x c ,y c ) in the Y direction;

[0106] Ang Ω represents the direction of the circular region Ω centered on the feature point (x c ,y c );

[0107] ξ represents a very small positive real number

[0108] The amplitude Mag c and direction Ang c of the circular region Ω where the feature point (x Ω , y Ω ) is located are obtained from equations (9) to (10). The range of the direction Ang Ω is set to [0, 360°), and it is divided into 36 intervals at intervals of 10°. The amplitude features and direction features of each interval are statistically analyzed to obtain the direction histogram distribution; the peak direction of the direction histogram is selected as the main direction of the feature point (x c , y c ).

[0109] A feature descriptor of the feature point (x c , y c ) is constructed using GLOH-like logarithmic polar coordinates: the radial direction of the circular region Ω where the feature point (x c , y c ) is located is divided into n r parts, and the angular direction is divided into n θ parts, resulting in (n r - 1)×n θ + 1 logarithmic polar coordinate sub-regions. The n bin direction amplitude features and angular features of each logarithmic polar coordinate sub-region are statistically analyzed to obtain the feature descriptor of the feature point (x c , y c ). By analogy, the feature descriptors corresponding to the remaining feature points on the simulated affine image S' are obtained.

[0110] The dimension of the feature descriptor corresponding to each feature point is n bin ×((n r - 1)×nθ + 1).

[0111] Then, the feature descriptors of all feature points on the simulated affine image S' are combined to obtain the feature descriptor Des S' of the simulated affine image S'.

[0112] Then, referring to the method for obtaining the feature descriptor Des S' of the simulated affine image S', the feature descriptor Des R of the orthoimage R is obtained.

[0113] Step 2.4, Feature Registration and Error Elimination

[0114] First, the feature descriptor Des S ' of the simulated affine image S′ and the feature descriptor Des RThe nearest neighbor distance ratio algorithm is used for initial registration, and the random sample consensus algorithm is used for fast outlier rejection to obtain the initial registration corresponding points of the simulated affine image S' and the orthophoto R; then the fast consensus sampling method is used to reject the incorrect registrations in the initial registration. When the residual of the corresponding points is less than 3 pixels, they are the homologous points of the correct registration of the simulated affine image S' and the orthophoto R.

[0115] Step 3. Image correction

[0116] Let the correctly registered feature point (x p , y p ) in the simulated affine image S' and the correctly registered feature point (x', p , y' p ) in the orthophoto R be a pair of homologous points. Rotate the correctly registered feature point (x p , y p ) in the simulated affine image S' according to the corresponding tilt T and the rotation angle in the longitude direction to find the corresponding feature point in the large tilt angle image S The homologous points of the large tilt angle image S and the orthophoto R are [[ID=, y' p , y' p ). And so on, all the homologous points between the large tilt angle image S and the orthophoto R are obtained.

[0117] Then, the affine transformation parameters Ψ between the large tilt angle image S and the orthophoto R are constructed through all the homologous points, and the large tilt angle image S and the orthophoto R are corrected according to the affine transformation parameters Ψ to complete the registration process.

[0118] Due to the adoption of the above technical solution, the present invention has the following positive effects compared with the prior art:

[0119] In step 1 of the present invention, the ASIFT algorithm is used to simulate the affine transformation perspective of the large tilt angle aerial image, which can better adapt to the large angle rotation problem. In step 2.1 of the present invention, the nonlinear diffusion filter is used to construct the nonlinear scale space. Compared with the traditional Gaussian pyramid scale space, the nonlinear scale space can retain the image boundary texture features while smoothing the image noise, and better cope with the problems of ground object geometric deformation, resolution difference, and ground object mutual occlusion existing in the large tilt angle image.

[0120] In step 2.2 of the present invention, PC is used to replace the traditional gradient information, and PC is combined with the Sobel operator to construct a feature point metric matrix to address the problem of complex ground object information in aerial images, and to extract richer and more stable structural features and point features of aerial images; in step 2.3, PC and the Sobel operator are used to calculate the amplitude and direction of feature points, and then a GLOH-like feature descriptor operator is constructed. Compared with the traditional SIFT-like feature descriptor operator, the amplitude and direction based on PC are invariant to image noise and brightness changes compared to those based on gradients, and can better represent the structural features and direction information of the image; in addition, the circular structure GLOH-like feature descriptor operator is more robust and stable than the square structure SIFT-like feature descriptor operator, resulting in a small registration error, a high registration success rate, and a stable registration result.

[0121] Therefore, the present invention has the characteristics of small registration error, good anti-noise performance, high success rate, and stable result. BRIEF DESCRIPTION OF THE DRAWINGS

[0122] Figure 1 It is a simulation diagram of longitude and latitude parameter sampling of the present invention;

[0123] Figure 2 It is a large tilt aerial image S of the present invention;

[0124] Figure 3 is Figure 2 the simulated affine image of the large tilt aerial image S shown;

[0125] Figure 4 It is a sub-simulated affine image of the present invention of the base space

[0126] Figure 5 It is a sub-simulated affine image of the present invention of the 3 sub-spaces;

[0127] Figure 6 It is a sub-simulated affine image of the present invention of the non-linear scale space;

[0128] Figure 7 It is the non-linear scale space of an orthoimage R of the present invention;

[0129] Figure 8 is Figure 6 the sub-simulated affine image shown of the PC feature map corresponding to the non-linear scale space;

[0130] Figure 9 is Figure 6 the sub-simulated affine image shown The feature point metric matrix M of the non - linear scale space i ;

[0131] Figure 10 is the process diagram for constructing a feature description operator of the present invention;

[0132] Figure 11 are the corresponding points of the correctly registered simulated affine image S' and ortho - image R of the present invention;

[0133] Figure 12 are the correctly registered feature points in the simulated affine image S' of the present invention and the corresponding feature points in the large - tilt image S;

[0134] Figure 13 is Figure 2 all the corresponding points between the large - tilt image S and the ortho - image R shown;

[0135] Figure 14 is Figure 2 the registration correction result of the large - tilt image S and the ortho - image R shown. Specific implementation manner

[0136] The present invention will be further described below in conjunction with the accompanying drawings and specific implementation manners, which is not a limitation of its protection scope.

[0137] Example 1

[0138] A method for registering large - tilt aerial images and ortho - images. The steps of the method in this embodiment are as follows:

[0139] Step 1: Use the ASIFT algorithm to construct the simulated affine image S' of the large - tilt aerial image S

[0140] Step 1.1: Parameter sampling

[0141] As Figure 1 shown, according to the ASIFT algorithm, in a hemispherical space, the viewing angle of the large - tilt aerial image S is simulated by the direction parameters of the camera in the longitude - direction angle φ and the latitude - direction angle θ; the longitude - direction angle φ and the latitude - direction angle θ are discretized respectively. Among them:

[0142] The latitude - direction angle θ adopts the geometric - series discretization method:

[0143] T n = a n (1) In formula (1):

[0144] T n represents the tilt of the camera at the nth sampling series in the latitude direction;

[0145] n represents the sampling level of the camera in the latitude direction, n = 1,..., N;

[0146] a represents the sampling interval of the camera in the latitude direction, a > 1.

[0147] The longitude direction angle φ adopts an arithmetic progression discretization method:

[0148]

[0149] In formula (2):

[0150] represents the inclination T of the camera in the latitude direction n under which the rotation angle at the k-th sampling level in the longitude direction;

[0151] b / T n represents the sampling interval in the longitude direction under the inclination T of the camera in the latitude direction, n where b is an angle constant;

[0152] b is an angle constant;

[0153] k represents the sampling level in the longitude direction under the inclination T of the camera in the latitude direction n where k = 1,..., K,

[0154] In this embodiment: N = 2; n = 1, 2; a = 1.2; b = 72°;

[0155] When n = 1, T1 = a = 1.2,

[0156] where, when n = 1, k = 1,

[0157] when n = 1, k = 2,

[0158] When n = 2, T2 = a 2 = 1.44,

[0159] where, when n = 2, k = 1,

[0160] when n = 2, k = 2,

[0161] when n = 2, k = 3,

[0162] Step 1.2, construct the simulated affine image S' of the large tilt angle aerial image S

[0163] Assume the camera's tilt in the latitude direction is T n , the rotation angle in the longitude direction is Simulated affine image of high-angle aerial image S The construction process is as follows: first, perform a bilinear interpolation rotation operation on the large-angle aerial image S to obtain a rotated image; then use an anti-aliasing filter to convolve the rotated image to obtain a filtered rotated image; then perform a tilt operation on the filtered rotated image to obtain the tilt of the camera in the latitude direction, which is T. n and the rotation angle in the longitude direction is Simulated affine image of high-angle aerial image S

[0164]

[0165] In formula (3):

[0166] Indicates that the camera's tilt in the latitude direction is T n The rotation angle in the longitude direction is The simulated affine image of the high-angle aerial image S at time ;

[0167] Tilt means to perform a tilt operation on the rotated image after filtering by the anti-aliasing filter;

[0168] G δ Represents the anti-aliasing filter, where δ represents the kernel standard deviation of the anti-aliasing filter, δ = 0.8T n ;

[0169] * indicates convolution operation;

[0170] Rot represents the bilinear interpolation rotation operation on the high-angle aerial image S;

[0171] S represents high-angle aerial imagery;

[0172] Indicates the tilt of the camera in the latitude direction T n The rotation angle of the k-th sampling level in the longitude direction;

[0173] T n Indicates the tilt of the camera at the nth sampling level in the latitude direction.

[0174] By analogy, the rest of the simulated affine images of the high-angle aerial image S are obtained. Then the simulated affine image of the high-angle aerial image S is

[0175] For the sake of simplicity of description: hereinafter, the "simulated affine image S' of the large tilt angle aerial image S" will be simply referred to as "simulated affine image S'"; and the "tilt angle of the camera in the latitude direction is T n and the rotation angle in the longitude direction is when the simulated affine image of the large tilt angle aerial image S " will be simply referred to as "sub-simulated affine image ".

[0176] In this embodiment:

[0177] The large tilt angle aerial image S is as Figure 2 shown. The large tilt angle aerial image S has a large tilt and rotation angle in the horizontal direction.

[0178] The simulated affine image of the large tilt angle aerial image S is as Figure 3 shown. Figure 3 Among them, from left to right are: the sub-simulated affine image with a tilt angle of T1 and a rotation angle of the sub-simulated affine image with a tilt angle of T1 and a rotation angle of the sub-simulated affine image with a tilt angle of T1 and a rotation angle of the sub-simulated affine image with a tilt angle of T2 and a rotation angle of the sub-simulated affine image with a tilt angle of T2 and a rotation angle of the sub-simulated affine image with a tilt angle of T2 and a rotation angle of the sub-simulated affine image with a tilt angle of T2 and a rotation angle of the sub-simulated affine image with a tilt angle of T2 and a rotation angle of the sub-simulated affine image with a tilt angle of T2 and a rotation angle of the sub-simulated affine image with a tilt angle of T2 and a rotation angle of

[0179] Step 2: Register the simulated affine image S' and the orthoimage R

[0180] Step 2.1: Construct a non-linear scale space

[0181] Discretize the non-linear scale space of the sub-simulated affine image into 1 base space and L sub-spaces. After discretization, the scale corresponding to the non-linear scale space of the sub-simulated affine image is:

[0182]

[0183] In formula (4):

[0184] i represents the i-th layer of the non-linear scale space of the sub-simulated affine image , i = 0, 1, 2,...., L;

[0185] In this embodiment: L = 3, i = 0, 1, 2, 3;

[0186] σ i represents the scale of the i-th layer of the non-linear scale space of the sub-simulation affine image. In this embodiment: When i = 0,

[0187] When i = 1,

[0188] When i = 2,

[0189] When i = 3,

[0190] When i = 3,

[0191] σ0 represents the initial scale of the non-linear scale space of the sub-simulation affine image The initial scale σ0 is a constant. In this embodiment: σ0 = 1.6;

[0192] L represents the number of sub-spaces of the sub-simulation affine image In this embodiment: L = 3.

[0193] Let the non-linear scale space of the sub-simulation affine image be: Where: represents the base space of the sub-simulation affine image represents the base space of the sub-simulation affine image successively represents the first sub-space, the second sub-space and the third sub-space of the sub-simulation affine image respectively.

[0194] The construction process of the base space of the sub-simulation affine image is as follows:

[0195] Use a Gaussian filter to perform smoothing filtering on the sub-simulation affine image The standard deviation of the Gaussian kernel of the Gaussian filter is the same as the initial scale σ0 of the non-linear scale space of the sub-simulation affine image In this embodiment: σ0 = 1.6. Obtain the base space Figure 4 of the sub-simulation affine image as shown

[0196] The construction process of the j-th sub-space of the sub-simulation affine image is as follows:

[0197] First, downsample the j-1-th layer of the non-linear scale space of the sub-simulation affine image at a sampling rate of 1 / σ0 = 1 / 1.6 = 0.625 to obtain a sampling space Then, use anisotropic diffusion filtering on the sampling space to perform filtering and obtain the j-th subspace of the sub-simulation affine image

[0198]

[0199] In Equation (5):

[0200] denotes the anisotropic diffusion filtering space obtained under the time step condition;

[0201] j represents the j-th of the subspaces of the sub-simulation affine image , where j = 1, 2,..., L;

[0202] In this embodiment: L = 3, j = 1, 2, 3;

[0203] ↓ represents the downsampling operation;

[0204] denotes downsampling the (j - 1)-th layer of the non-linear scale space of the sub-simulation affine image to obtain the sampling space of

[0205] ;

[0206] t represents the time measure;

[0207] div represents the divergence operator;

[0208] D(x, y, t) represents the diffusion coefficient, and the diffusion coefficient depends on the gradient norm

[0209] denotes the gradient operator;

[0210] denotes the j-th subspace of the sub-simulation affine image ;

[0211] I d denotes the identity matrix;

[0212] σ j denotes the scale of the j-th layer of the non-linear scale space of the sub-simulation affine image. In this embodiment:

[0213] When j = 1,

[0214] When j = 2,

[0215] ​​When j = 3,

[0216] σ j-1 represents the scale of the (j - 1)-th layer non-linear scale space of the sub-simulated affine image. In this embodiment: When j = 1,

[0217] When j = 1,

[0218] When j = 2,

[0219] When j = 3,

[0220] l represents the direction axis of anisotropic diffusion;

[0221] m represents the number of direction axes of anisotropic diffusion. In this embodiment: m = 2;

[0222] A l represents the diffusion coefficient matrix of the sampling space along the direction axis l of anisotropic diffusion. A l is the discrete form of the diffusion coefficient D.

[0223] According to Equation (5), L = 3, j = 1, 2, 3; the three sub-spaces of the sub-simulated affine image Figure 5 as shown in are obtained. From left to right in Figure 5 they represent the first sub-space the second sub-space and the third sub-space of the sub-simulated affine image Then, as shown in Figure 6 the non-linear scale space of the sub-simulated affine image is: And so on, the non-linear scale spaces of the remaining sub-simulated affine images in the simulated affine image S' are obtained.

[0224] By referring to the method of the non-linear scale spaces of the remaining sub-simulated affine images in the simulated affine image S', the non-linear scale space of the orthographic image R as shown in Figure 7 is: R0, R1,..., R L .

[0225] Step 2.2, Feature Point Detection

[0226] For the sake of convenience of description, the phase consistency is represented by "PC".

[0227] Feature point detection is performed on the sub-simulated affine image The feature point detection process is as follows:

[0228] First, use PC to extract the structural features of the i-th layer non-linear scale space image of the sub-simulated affine image to obtain the sub-simulated affine image of the i-th layer non-linear scale space PC feature map; then use the Sobel operator templates in the X and Y directions to convolve the sub-simulated affine image of the i-th layer non-linear scale space PC feature map to obtain the first-order derivatives of the sub-simulated affine image of the i-th layer non-linear scale space PC feature map in the X and Y directions:

[0229]

[0230]

[0231] In equations (6) to (7):

[0232] i represents the i-th layer of the non-linear scale space of the sub-simulated affine image where i = 0, 1, 2,..., L; in this example: L = 3, i = 0, 1, 2, 3;

[0233] represents the first-order derivative of the sub-simulated affine image of the i-th layer non-linear scale space PC feature map in the X direction;

[0234] represents the first-order derivative of the sub-simulated affine image of the i-th layer non-linear scale space PC feature map in the Y direction;

[0235] PC i represents the sub-simulated affine image of the i-th layer non-linear scale space PC feature map;

[0236] * represents the convolution operation;

[0237] Γ x represents the Sobel operator template in the X direction;

[0238] Γ y represents the Sobel operator template in the Y direction.

[0239] Then, use the first-order derivatives of the sub-simulated affine image of the i-th layer non-linear scale space PC feature map in the X and Y directions to construct the feature point metric matrix M of the i-th layer non-linear scale space of the sub-simulated affine image i :

[0240] ​

[0241] In formula (8):

[0242] i represents the i-th layer of the non-linear scale space of the sub-simulated affine image , where i = 0, 1, 2,..., L; in this embodiment: L = 3, i = 0, 1, 2, 3;

[0243] M i represents the feature point metric matrix of the i-th layer of the non-linear scale space of the sub-simulated affine image ;

[0244] represents the Gaussian kernel, and the standard deviation of the Gaussian kernel is where σ i represents the scale of the i-th layer of the scale space of the sub-simulated affine image ; in this embodiment:

[0245] When i = 0,

[0246] When i = 1,

[0247] When i = 2,

[0248] When i = 3,

[0249] * represents the convolution operation;

[0250] represents the first derivative of the PC feature map of the i-th layer of the non-linear scale space of the sub-simulated affine image in the X direction;

[0251] represents the first derivative of the PC feature map of the i-th layer of the non-linear scale space of the sub-simulated affine image in the Y direction.

[0252] Finally, perform eigen-decomposition on the feature point metric matrix M of the i-th layer of the non-linear scale space of the sub-simulated affine image to obtain the first eigenvalue λ1 and the second eigenvalue λ2 of the feature point metric matrix M i of the i-th layer of the non-linear scale space of the sub-simulated affine image. The feature point intensity map of the i-th layer of the non-linear scale space of the sub-simulated affine image is: R i = min(λ1, λ2), and the points satisfying R i > r are the points of the sub-simulated affine image i ; i >r of the sub-simulated affine image The feature points of the i-th layer of the non-linear scale space, where r is a threshold constant. In this embodiment: r = 0.004.

[0253] In this embodiment:

[0254] Sub-simulated affine image The PC feature map corresponding to the non-linear scale space of Figure 8 is shown as Figure 8 In it, from left to right are: the PC feature map of the 0-th layer of the non-linear scale space of the sub-simulated affine image , the PC feature map of the 1-st layer of the non-linear scale space, the PC feature map of the 2-nd layer of the non-linear scale space, and the PC feature map of the 3-rd layer of the non-linear scale space.

[0255] Sub-simulated affine image The feature point metric matrix M of the non-linear scale space of i is shown as Figure 9 in Figure 9 In it, from left to right are: the feature point metric matrix M0 of the 0-th layer of the non-linear scale space of the sub-simulated affine image , the feature point metric matrix M1 of the 1-st layer of the non-linear scale space, the feature point metric matrix M2 of the 2-nd layer of the non-linear scale space, and the feature point metric matrix M3 of the 3-rd layer of the non-linear scale space.

[0256] [[ID=no]]In the feature points of each layer of the non-linear scale space of the obtained sub-simulated affine image , if there are more than 2 feature points at the same position in different layers of the non-linear scale space, only one of these feature points is retained; if there are feature points at different positions in different layers of the non-linear scale space, all of these feature points are retained, obtaining the feature point set of the sub-simulated affine image

[0257] And so on, perform feature point detection on the remaining sub-simulated affine images in the simulated affine image S'.

[0258] Combine the feature point sets of all sub-simulated affine images to obtain the feature point set

[0259] Refer to the method of the feature point set of the simulated affine image S' to obtain the feature point set Key R = {(x1', y1'),..., (x' m , y' m )} of the orthophoto R.

[0260] Step 2.3, Construction of feature description operator

[0261] For any feature point (x c , y c ) in the simulated affine image S', extract the amplitude Mag c and direction Ang c of the circular region Ω centered on the feature point (x Ω and direction Ang Ω :

[0262]

[0263]

[0264] In equations (9) to (10):

[0265] Ω represents a circular region centered on the feature point (x c , y c ), and the radius of the circular region Ω is d r . In this embodiment: d r = 54;

[0266] Mag Ω represents the amplitude of the circular region Ω centered on the feature point (x c , y c );

[0267] represents the first derivative of the PC feature map of the circular region Ω centered on the feature point (x c , y c ) in the X direction;

[0268] represents the first derivative of the PC feature map of the circular region Ω centered on the feature point (x c , y c ) in the Y direction;

[0269] Ang Ω represents the direction of the circular region Ω centered on the feature point (x c , y c );

[0270] ξ represents a very small positive real number. In this embodiment: ξ = 0.001.

[0271] From equations (9) to (10), obtain the amplitude Mag c and direction Ang c of the circular region Ω where the feature point (x Ω and direction Ang Ω is located. Let the direction Ang ΩThe range is set to [0, 360°), and then it is divided into 36 intervals at an interval of every 10°. The amplitude characteristics and direction characteristics of each interval are counted to obtain the direction histogram distribution; the peak direction of the direction histogram is selected as the main direction of the feature point (x c , y c ).

[0272] Construct the feature descriptor of the feature point (x c , y c ) using GLOH-like logarithmic polar coordinates as Figure 10 shown: Divide the polar radius direction of the circular region Ω where the feature point (x c , y c ) is located into n r = 5 parts, and divide the polar angle direction into n θ = 9 parts, obtaining (n r - 1) × n θ + 1 = 37 logarithmic polar coordinate sub-regions. Count the n bin = 8 direction amplitude characteristics and angle characteristics of each logarithmic polar coordinate sub-region to obtain the feature descriptor of the feature point (x c , y c ). By analogy, obtain the feature descriptors corresponding to the remaining feature points on the simulated affine image S'.

[0273] The dimension of each of the feature descriptors is n bin × ((n r - 1) × n θ + 1) = 296.

[0274] Then combine the feature descriptors of all feature points on the simulated affine image S' to obtain the feature descriptor Des S' of the simulated affine image S'.

[0275] Then, referring to the method for obtaining the feature descriptor Des S' of the simulated affine image S', obtain the feature descriptor Des R of the orthoimage R.

[0276] Step 2.4, Feature registration and error elimination

[0277] First, perform initial registration on the feature descriptor Des S' of the simulated affine image S' and the feature descriptor Des R of the orthoimage R using the nearest neighbor distance ratio algorithm, and perform fast outlier elimination using the random sample consensus algorithm to obtain the initial registration corresponding points of the simulated affine image S' and the orthoimage R; then use the fast consensus sampling method to eliminate the misregistrations in the initial registration. When the residual of the corresponding points is less than 3 pixels, such asFigure 11 As shown, they are the corresponding points of the simulated affine image S' and the orthoimage R that are correctly registered.

[0278] Step 3: Image correction

[0279] Let the correctly registered feature point (x p , y p ) in the simulated affine image S' and the correctly registered feature point (x' p , y' p ) in the orthoimage R be a pair of corresponding points. Rotate the correctly registered feature point (x p , y p ) in the simulated affine image S' according to the corresponding tilt angle T and the rotation angle in the longitude direction to find the corresponding feature point in the large tilt angle image S as shown in Figure to obtain the corresponding points of the large tilt angle image S and the orthoimage R as and (x' , y' p , y' p ). By analogy, all the corresponding points between the large tilt angle image S and the orthoimage R as shown in ​ are obtained.

[0280] Then, construct the affine transformation parameters Ψ between the large tilt angle image S and the orthoimage R through all the corresponding points, and correct the large tilt angle image S and the orthoimage R according to the affine transformation parameters Ψ. As shown in ​ , the registration process is completed. As can be seen from ​ , the large tilt angle image S and the orthoimage R can overlap well in the common feature area and fit perfectly at the area boundary, indicating that the registration error of this example is small, the anti-noise performance is good, the success rate is high, and the result is stable.

[0281] In this specific implementation manner, the ASIFT algorithm is adopted in step 1 to simulate the affine transformation perspective of the large tilt angle aerial image, which can better adapt to the large angle rotation problem. In step 2.1 of this specific implementation manner, the nonlinear diffusion filter is adopted to construct the nonlinear scale space. Compared with the traditional Gaussian pyramid scale space, the nonlinear scale space can retain the image boundary texture features while smoothing the image noise, and can better handle the problems of ground object geometric deformation, resolution difference, and ground object mutual occlusion existing in the large tilt angle image.

[0282] In this specific embodiment, in step 2.2, PC is used to replace the traditional gradient information, and PC is combined with the Sobel operator to construct a feature point metric matrix to address the problem of complex ground object information in aerial images and extract richer and more stable structural features and point features of aerial images; in step 2.3, PC and the Sobel operator are used to calculate the amplitude and direction of feature points, and then a GLOH-like feature descriptor operator is constructed. Compared with the traditional SIFT-like feature descriptor operator, the amplitude and direction based on PC are invariant to image noise and brightness changes compared to those based on gradients, and can better represent the structural features and direction information of images; in addition, the circular structure GLOH-like feature descriptor operator is more robust and stable than the square structure SIFT-like feature descriptor operator, resulting in a small registration error, a high registration success rate, and a stable registration result.

[0283] Therefore, this specific embodiment has the characteristics of small registration error, good anti-noise performance, high success rate, and stable results.

Claims

1. A method for registering high-angle aerial images with orthophotos, characterized in that The specific steps of the high-angle aerial image and orthophoto registration method are as follows: Step 1: Use the ASIFT algorithm to construct the simulated affine image S' of the high-angle aerial image S Step 1.1: Parameter sampling According to the ASIFT algorithm, in a hemispherical space, the viewing angle of the high-angle aerial image S is simulated by the directional parameters of the camera in the longitude direction angle φ and the latitude direction angle θ. The longitude direction angle φ and the latitude direction angle θ are discretized respectively, where: The latitude angle θ is discretized using geometric series: T n =a n (1) In formula (1): T n Indicates the tilt of the camera at the nth sampling level in the latitude direction, n represents the sampling level of the camera in the latitude direction, n=1,...,N, a represents the sampling interval of the camera in the latitude direction, a>1; The longitude angle φ is discretized using the arithmetic level: In formula (2): Indicates the tilt of the camera in the latitude direction T n The rotation angle of the k-th sampling level in the longitude direction is, b / T n Indicates the tilt of the camera in the latitude direction T n In the sampling interval in the longitude direction, b is an angle constant, k represents the tilt of the camera in the latitude direction T n Under the condition of , the number of sampling levels in the longitude direction is k=1,...,K, Step 1.2: Construct the simulated affine image S' of the high-angle aerial image S Assume the camera's tilt in the latitude direction is T n , the rotation angle in the longitude direction is Simulated affine image of high-angle aerial image S The construction process is as follows: first, perform a bilinear interpolation rotation operation on the large-angle aerial image S to obtain a rotated image; then use an anti-aliasing filter to convolve the rotated image to obtain a filtered rotated image; then perform a tilt operation on the filtered rotated image to obtain the tilt of the camera in the latitude direction, which is T. n and the rotation angle in the longitude direction is Simulated affine image of high-angle aerial image S In formula (3): Indicates that the camera's tilt in the latitude direction is T n The rotation angle in the longitude direction is The simulated affine image of the high-angle aerial image S at time , Tilt means to perform a tilt operation on the rotated image after filtering by the anti-aliasing filter. G δ Represents the anti-aliasing filter, where δ represents the kernel standard deviation of the anti-aliasing filter, δ = 0.8T n , * indicates the convolution operation, Rot represents the bilinear interpolation rotation operation on the high-angle aerial image S. S represents high-angle aerial imagery; By analogy, the rest of the simulated affine images of the high-angle aerial image S are obtained. Then the simulated affine image of the high-angle aerial image S is For the sake of simplicity, the following will refer to the simulated affine image S' of the large-angle aerial image S as "simulated affine image S'"; the camera's tilt in the latitude direction is T n and the rotation angle in the longitude direction is Simulated affine image of high-angle aerial image S "Sub-simulated affine image" Step 2: Register the simulated affine image S' and the orthophoto image R Step 2.1: Nonlinear scale space construction Simulate affine images The nonlinear scale space is discretized into 1 base space and L sub-level spaces. After discretization, the sub-simulated affine image The scale corresponding to the nonlinear scale space is: In formula (4): i represents the sub-simulated affine image The i-th layer of the nonlinear scale space, i=0,1,2,....,L, σ i Represents a sub-simulated affine image The scale of the i-th layer nonlinear scale space, σ0 represents the sub-simulated affine image The initial scale of the nonlinear scale space, σ0 is a constant, L represents the sub-simulated affine image The number of sub-spaces; Assume that the sub-simulation affine image The nonlinear scale space is: in Represents a sub-simulated affine image grassroots space, Sub-simulated affine images are represented in sequence The 1st sub-level space, ..., the Lth sub-level space; Sub-simulated affine image Grassroots space The build process is as follows: Simulating affine images using Gaussian filter pairs Perform smoothing filtering, Gaussian kernel standard deviation of Gaussian filter and sub-simulated affine image The initial scale σ0 of the nonlinear scale space is the same, and the sub-simulated affine image is obtained Grassroots space Sub-simulated affine image The j-th sub-space The build process is as follows: First, simulate the affine image with a sampling rate of 1 / σ0 The j-1th layer nonlinear scale space is downsampled to obtain the sampling space Then use anisotropic diffusion filtering to filter the sampling space Filter and obtain the sub-simulated affine image The j-th sub-space In formula (5): Indicates the time step The anisotropic diffusion filter space obtained under the conditions, j represents the sub-simulated affine image The jth sub-space of , j = 1, 2, ..., L, ↓ represents the downsampling operation, Represents a pair of simulated affine images The sampling space obtained by downsampling the j-1th layer nonlinear scale space is t represents the time measurement, div represents the divergence operator, D(x,y,t) represents the diffusion coefficient, which depends on the gradient norm represents the gradient operator, Represents a sub-simulated affine image The j-th sub-level space of I d represents the identity matrix, σ j Represents a sub-simulated affine image The scale of the j-th layer nonlinear scale space, σ j-1 Represents a sub-simulated affine image The scale of the j-1th layer nonlinear scale space, l represents the direction axis of anisotropic diffusion, m represents the number of directional axes of anisotropic diffusion, A l Represents the diffusion coefficient matrix of the sampling space along the direction axis l of anisotropic diffusion, A l is the discrete form of the diffusion coefficient D; According to formula (5), the sub-simulated affine image is obtained L sub-spaces Zi simulated affine image The nonlinear scale space is: By analogy, the nonlinear scale space of the remaining sub-simulated affine images in the simulated affine image S' is obtained; Referring to the method of obtaining the nonlinear scale space of the other sub-simulated affine images in the simulated affine image S', the nonlinear scale space of the orthophoto R is obtained as follows: R0, R1, ..., R L ; Step 2.2, feature point detection For the convenience of description, phase consistency is represented by "PC"; Pair simulated affine image Perform feature point detection, the feature point detection process is as follows: First use PC to extract the simulated affine image The structural characteristics of the i-th layer nonlinear scale space image are obtained by sub-simulating the affine image The PC feature map of the i-th layer nonlinear scale space; then use the Sobel operator template pair in the X direction and the Y direction to simulate the affine image Convolution is performed on the PC feature map of the i-th layer nonlinear scale space to obtain the sub-simulated affine image The first-order derivatives of the PC feature map of the i-th layer nonlinear scale space in the X and Y directions are: In formulas (6) to (7): i represents the sub-simulated affine image The i-th layer of the nonlinear scale space, i=0,1,2,...,L, Represents a sub-simulated affine image The first-order derivative of the PC feature map of the i-th layer nonlinear scale space in the X direction, Represents a sub-simulated affine image The first-order derivative of the PC feature map of the i-th layer nonlinear scale space in the Y direction, PC i Represents a sub-simulated affine image The PC feature map of the i-th layer nonlinear scale space, * indicates the convolution operation, Γ x Represents the Sobel operator template in the X direction, Γ y Represents the Sobel operator template in the Y direction; Then use the obtained sub-simulation affine image The first-order derivatives of the PC feature map of the i-th layer nonlinear scale space in the X and Y directions are used to construct the sub-simulated affine image The feature point measurement matrix M of the i-th layer nonlinear scale space i : In formula (8): i represents the sub-simulated affine image The i-th layer of the nonlinear scale space, i=0,1,2,....,L, M i Represents a sub-simulated affine image The feature point measurement matrix of the i-th layer nonlinear scale space, represents the Gaussian kernel, and the standard deviation of the Gaussian kernel is where σ i Represents a sub-simulated affine image The scale of the i-th level scale space, * indicates convolution operation; Finally, the pair simulates an affine image The feature point measurement matrix M of the i-th layer nonlinear scale space i Perform feature decomposition to obtain the sub-simulated affine image The feature point measurement matrix M of the i-th layer nonlinear scale space i The first eigenvalue λ1 and the second eigenvalue λ2 of the sub-simulated affine image The intensity map of the feature point in the i-th layer nonlinear scale space is: R i =min(λ1,λ2), satisfying R i >r is a sub-simulated affine image The feature points of the i-th layer nonlinear scale space, r is a threshold constant; In the obtained sub-simulated affine image Among the feature points of each layer of nonlinear scale space, if there are more than two feature points at the same position in the nonlinear scale space of different layers, only one of these feature points is retained; if there are feature points at different positions in the nonlinear scale space of different layers, all of these feature points are retained to obtain a sub-simulated affine image. A set of feature points Similarly, feature point detection is performed on the remaining sub-simulated affine images in the simulated affine image S'; Combine the feature point sets of all sub-simulated affine images to obtain the feature point set of the simulated affine image S' Refer to the method of simulating the feature point set of the affine image S' to obtain the feature point set Key of the orthophoto R R ={(x1',y1'),...,(x' m ,y' m )}; Step 2.3: Feature description operator construction For any feature point (x c ,y c ), extract the feature points (x c ,y c ) is the amplitude of the circular area Ω centered at Ω and direction Ang Ω : In formulas (9) and (10): Mag Ω Indicates that the feature point (x c ,y c ) is the amplitude of the circular area Ω centered at Ω represents the characteristic point (x c ,y c ) is the center of the circular area, the radius of the circular area Ω is d r , Indicates that the feature point (x c ,y c ) is the first-order derivative of the PC feature map in the X direction of the circular area Ω centered at Indicates that the feature point (x c ,y c ) is the first derivative of the PC feature map in the Y direction of the circular area Ω centered at Ang Ω Indicates that the feature point (x c ,y c ) is the direction of the circular area Ω centered, ξ represents a very small positive real number; The characteristic point (x c ,y c ) in the circular area Ω where the amplitude Mag Ω and direction Ang Ω , change the direction to Ang Ω The range is set to [0,360°), and then divided into 36 intervals at intervals of 10°; the amplitude characteristics and direction characteristics of each interval are counted to obtain the direction histogram distribution; the peak direction of the direction histogram is selected as the feature point (x c ,y c )’s main direction; The feature points (x c ,y c ) feature description operator: the feature point (x c ,y c ) is divided into the polar direction of the circular area Ω r The polar angle direction is divided into n θ , get (n r -1)×n θ+ 1 logarithmic polar coordinate sub-region; count the n of each logarithmic polar coordinate sub-region bin Direction amplitude features and angle features are obtained to obtain the feature point (x c ,y c ) feature description operator; and so on, the feature description operators corresponding to the remaining feature points on the simulated affine image S' are obtained; The dimension of the feature description operator corresponding to each feature point is n bin ×((n r -1)×n θ +1); Then combine the feature description operators of all feature points on the simulated affine image S' to obtain the feature description operator Des of the simulated affine image S' S' ; Then refer to the feature description operator Des of the simulated affine image S' S' The method is used to obtain the feature description operator Des of the orthophoto R. R ; Step 2.4: Feature registration and error elimination First, the feature description operator Des of the simulated affine image S' S' and the feature description operator Des of the orthophoto R R The nearest neighbor distance ratio algorithm is used for initial registration, and the random sampling consistency algorithm is used for fast outlier removal to obtain the initial registration corresponding points of the simulated affine image S' and the orthophoto image R. The fast consistency sampling method is then used to remove the incorrect registration in the initial registration. When the residual of the corresponding point is less than 3 pixels, it is the same name point of the correct registration of the simulated affine image S' and the orthophoto image R. Step 3: Image Correction Suppose the correctly registered feature points (x p ,y p ) and the correctly registered feature points (x' p ,y' p ) are a pair of points with the same name, and the correctly registered feature points (x p ,y p ) Rotate according to the corresponding tilt T and longitude angle Find the feature points corresponding to the large-angle image S The same-name points of the large-angle image S and the orthophoto image R are obtained as and (x' p ,y' p ); By analogy, all the same-name points between the high-angle image S and the orthophoto image R are obtained; then, the affine transformation parameter Ψ between the high-angle image S and the orthophoto image R is constructed through all the same-name points, and the high-angle image S and the orthophoto image R are corrected according to the affine transformation parameter Ψ to complete the registration process.

Images