CN104766323B - A kind of Point matching method of remote sensing images - Google Patents

Abstract

The invention discloses a point matching method of a remote sensing image, which is used for screening out a correct point matching set from an initial matching point set of a remote sensing image with a large background change. This method is based on the similarity of the global structure formed by the correct matching points in the two images. First, a number of corresponding triangle pairs are composed of three randomly selected point pairs, and the difference between the corresponding angles and the corresponding side length ratios of the two triangles is calculated, and the spatial rotation order of the three vertices relative to the reference point is compared, and the similar triangles are selected. Conditional pair of triangles. Then count the occurrence frequency of each point pair in these similar triangles, build a histogram, set a threshold, and remove the point pairs whose frequency is less than a certain threshold. Finally, a spatial consistency measure is used to remove mismatched points whose coordinates are very similar. The invention improves the point matching precision of the remote sensing image, and at the same time improves the registration precision of the remote sensing image.

Description

A Point Matching Method of Remote Sensing Image

technical field

The invention relates to the technical field of remote sensing image processing, in particular to a point matching method of remote sensing images.

Background technique

Image registration is one of the most basic techniques in image processing, and it has been widely used in fields such as medicine, pattern recognition, remote sensing and computer vision. Registration is to align different images so that all the features can be expressed using the same coordinate system. Feature matching is the most important link in image registration. If there are too many mismatches in the feature matching set, wrong transformation parameters will be generated. Therefore a robust and accurate feature matching method is necessary. Most of the currently used feature matching methods are point matching methods. Point matching methods are mainly divided into two types: 1) matching using local gray information; 2) matching using local neighborhood structure or global structure. SIFT is a point feature extraction and point matching method based on local description operators, which find matching point pairs according to the similarity of point descriptors. However, when this method is applied to remote sensing images, especially images before and after disasters, the local gray value changes greatly, and it is very unreliable to rely solely on local gray information, and it is easy to generate more mismatching points. RANSAC is a classic point matching method, which uses the initial matching point set to estimate the transformation parameters, and removes the false matches in the matching point set. However, this method does not work well when the contrast ratio of mismatched points is large. GTM uses the spatial relationship in the field of feature points K to find the correct match from the initial matching point set, but this method cannot rule out the wrong match with the same neighborhood structure. In addition, correct matching points surrounded by a large number of false matching points may also be removed by mistake. The WGTM algorithm is an improved version of GTM, which uses the angular distance between the line segments connecting the feature points as the weight. But it can only remove some of the false matches with similar neighborhood structures, and it cannot distinguish the correct matching points with different neighborhood structures. The Bi-SOGC algorithm is a method based on GTM and is mainly used in the field of remote sensing images. It uses ordered line segments instead of unordered line segments, the spatial order of feature points is used to remove false matches with similar neighborhood structures, and a recovery strategy is used to retrieve correct point pairs that were removed by mistake. This method is effective for some remote sensing images, but it cannot deal with isolated correct matching points without any feature points in the neighborhood.

Some of the existing point matching methods are based on local gray information, and some are based on local neighborhood structure. However, when it is applied to remote sensing image registration, especially the remote sensing images before and after disasters, due to the large background changes, there are often many false matches in the initial matching set, and there are few correct matches. At this time, the distribution of correct matches is sparse, and the neighborhood can There is less correct structural information to use, and only using local grayscale information or neighborhood structural information often cannot achieve satisfactory results. The present invention (Triangle Transformation Matching Method TTM) utilizes the similarity of the global structure formed by correct matching points to remove false matches and retain correct matches. Although the gray information of the remote sensing images before and after the disaster changes greatly, the global structure composed of correct matching points is still similar. The invention uses the triangle as the primitive of the global structure to compare its similarity, and finally extracts the correct matching point set.

Contents of the invention

The purpose of the present invention is to provide a point matching method of remote sensing images, which improves the point matching accuracy of remote sensing images and simultaneously improves the registration accuracy of remote sensing images.

The technical solution adopted in the present invention is: a method for point matching of remote sensing images, the method comprising the following steps:

Step (1): The initial matching point set is two corresponding point sets extracted from the input in the reference image and the sensing image, in the initial matching point set P={p _i } and P'={p _i '} , i=1,2,3...N, N is the number of points in the point set, three point pairs are selected arbitrarily to form many corresponding triangle pairs, and the three point pairs selected arbitrarily are expressed as (p _i ,p _j ,p _k ) and (p _i ', p _j ', p _k ');

Step (2): Calculate the corresponding angle difference and the corresponding side length ratio difference between the two triangles formed:

Step (3): and Indicates the horizontal and vertical coordinates of points p _i and p _i ', respectively set the reference points p _d and p _d ' in the two triangles, and the coordinate calculation formula is as follows:

Step (4): Set triangle similarity threshold θ _b = 5°, l _b = 0.1, and then according to the order consistency and conditions of the three vertices relative to the reference point △θ ₁ , △θ ₂ , △θ ₃ ∈[0 ,θ _b ], △l ₁ ,△l ₂ ,△l ₃ ∈[0,l _b ] to filter out the qualified triangle pairs;

Step (5): In the selected triangle pair, count the number of occurrences of each point pair and establish a frequency histogram;

Step (6): Set the threshold k _p according to the experimental analysis, and remove the point pairs whose occurrence times are less than k _p ;

Step (7): Calculate the model conversion parameters according to all point pairs at this time, the point pairs are expressed as P _o and P _o ', and the number is n;

T=argminE(T(θ)) (10)

in, Represents the root mean square error, T(θ) represents the conversion parameters obtained by the point set P _o and P _o 'using the least square method, d(i)=||T(θ)(p _i ')-p _i ||, i=1,2,3,...n represents the conversion error of each point pair;

Step (8): Set the threshold E _t , when E(T(θ))>E _t , the point pair with the largest d(i) value will be removed, and the point set will be updated at the same time, when the condition E(T (θ))≤E _t , the iteration ends, and two correct matching point sets are obtained at this time.

Among them, the method is used to filter out the correct point matching set from the initial matching point set of remote sensing images with large background changes.

The advantage of the present invention compared with prior art is:

(1) The present invention utilizes the global structure as a comparison element and uses a triangle as a primitive, which is easy to calculate and can ensure the stability of the matching result.

(2) The present invention ensures the robustness of the results by establishing a histogram and performing screening according to the statistical results.

(3) The present invention uses the global consistency criterion to remove false matches with similar coordinates and improve the correct match rate.

Description of drawings

Figure 1 shows the global structure composed of correctly matched feature points, where (a) is the image before the disaster and the global structure of its feature points, and (b) is the image after the disaster and its corresponding correctly matched feature points Constituted global structure diagram;

Fig. 2 is the whole flowchart of the point matching method of a kind of remote sensing image of the present invention;

Figure 3 is an example of a mismatch point;

Figure 4 is a graph showing the variation of matching results with _kp ;

Figure 5 (a) is the initial matching result, (b) is the TTM matching result;

Fig. 6 is the frequency statistics histogram (being frequency histogram or frequency distribution histogram) of feature point;

Figure 7 shows (a) the initial matching result, (b) the RANSAC matching result, (c) the GTM matching result, (d) the Bi-SOGC matching result, and (e) the TTM matching result;

Fig. 8 is a graph of matching accuracy results of different spectral images;

Figure 9(a) is the initial matching result, (b) is the RANSAC matching result, (c) is the GTM matching result, (d) is the Bi-SOGC matching result, (e) is the TTM matching result;

Figure 10 shows the results of image matching accuracy before and after the disaster.

Detailed ways

The present invention will be further described below in conjunction with the accompanying drawings and specific embodiments.

The input of the method of the present invention is two corresponding point sets extracted in the reference image and the sensing image, which can be expressed as P={p _i } and P'={p _i '}(i=1,2, 3...N), (p _i corresponds to p _i '). In the initial set of matching points, there are many mismatching points, especially in remote sensing images with background changes. The TTM algorithm is based on the similarity of the overall structure composed of correct points, and it is used to remove most of the mismatch points and keep all the correct points. Since there are often affine transformations between different remote sensing images taken in the same area, triangles are chosen instead of quadrilaterals or other polygons as the basic element of the overall structure because of their better shape retention. Figure 1 shows two remote sensing images before and after the disaster, and it can be seen that the overall structure composed of correct matching points is similar. As shown in FIG. 2 , the method includes three steps: Step 1, feature extraction corresponding to the triangle; Step 2, building a frequency histogram; Step 3, removing residual mismatch points. Below in conjunction with specific accompanying drawing, the invention is described in detail as follows:

1. Feature extraction corresponding to the triangle

Three matching point pairs are randomly selected in the initial matching point set, which are described as (p _i , p _j , p _k ) and (p _i ', p _j ', p _k '). Then perform the following steps:

(1) Form two triangles in two images by (p _i , p _j , p _k ) and (p _i ', p _j ', p _k '), and then calculate the difference between the corresponding angle and the corresponding side length ratio :

If the three point pairs are all correct matches, the conditions △θ ₁ , △θ ₂ , △θ ₃ ∈ [0,θ _b ], △l ₁ , △l ₂ , △l ₃ ∈ [0,l _b ]. Three-point pair combinations that meet the conditions will be recorded.

(2) As shown in Figure 3, (p _a , p _a ') and (p _b , p _b ') are two correct matches, and (p _c , p _c ') is a wrong match. But this case still satisfies condition (1), and (p _c ,p _c ') will be incorrectly recorded as a correct match. Therefore, define the reference point (p _d ,p _d '), and the coordinate calculation formula is as follows:

In Figure 3, let the vector Rotate clockwise around the point p _d , p _d ', if the three feature point pairs are all correct matches, the order of the vector passing through the three vertices should be the same. As shown in Figure 3, the sequence of points passed by the left image is (p _a , p _c , p _b ), while the sequence of points passed by the right image is (p _a ', p _b ', p _c '), so (p _c ,p _c ') will not be considered a correct match. The combination of three point pairs satisfying both condition 1 and condition 2 will be selected into subset S, and points like (p _a ,p _a ') and (p _b ,p _b ') will be combined with another correct match Selected when. Due to the similarity of the structure formed by the correct matching points, the correct matching will be selected into S (as long as there are at least three pairs of correct matching), but some false matchings will also be wrongly selected, one with each point in the subset S The number-dependent histogram will be built in (b).

2. Create a frequency histogram

After step 1 and step 2 of part (1) are executed for each three-point pair combination, the number of occurrences of each point in S will be counted. Theoretically, as long as there is a mismatch point in the combination, this combination will not be selected into S. Therefore, it can be seen that in S, the number of correct matches must be much greater than the number of wrong matches. The frequency histogram is built as follows:

hist _P ={k(p _i ),i=1,2,...,M} (5)

k(p _i ) represents the number of combinations containing points (p _i , p _i ') in the subset. M represents the number of all unique points in S. From this we can conclude that the maximum number of entries is:

k _max ＝max(k(p ₁ ),k(p ₂ ),...,k(p _M )) (6)

In order to determine the correct threshold k _p , the recall and precision accuracy of 20 pairs of images under different thresholds are shown in Figure 4.

The calculation formulas of recall and precision are as follows:

recall=(number of correct matching points extracted)/(all correct matching points) (7)

precision=(extracted correct matching points)/(all extracted points) (8)

If the number of correct point pairs exists is m, the number of correct point pairs in S can reach much larger than the number of mismatches. Combined with Figure 4, so set the threshold:

At the same time, the point pairs whose number is less than k _p are eliminated. Fig. 6 is a statistical diagram of the frequency (ie frequency) of all the feature points in Fig. 5. After verification, the point pairs whose frequency exceeds the horizontal line in the figure are all correct points, and no correct point is below the horizontal line.

3. Remove residual mismatches

Most of the mismatches will be removed in the first two parts. But there are also points similar to ( _pe , _pe ') in Fig. 3 . (p _e , _pe ') is a mismatch, but because their coordinates are very similar, such mismatches cannot be removed in the previous process. Therefore the spatial consistency measure is used to judge whether a point pair is mismatched or not. After part A and part B, two matching point sets P _o and P _o ' with n points are obtained, and the conversion parameter T between images can be calculated by the following formula:

T=argminE(T(θ)) (10)

Indicates the root mean square error (RMSE), and T(θ) indicates the conversion parameters of the point sets P _o and P _o ' using the least square method. d(i)=||T(θ)(p _i ')-p _i ||, i=1,2,3,...n.

In order to remove false matches with similar coordinates and obtain optimal transformation parameters, a threshold E _t is set. During each iteration, when E(T(θ))>E _t , the point pair with the largest d(i) value will be removed, and the point set will be updated at the same time. When the condition E(T(θ)) _≤Et is satisfied, the iteration ends. At this point, two correct matching point sets can be obtained.

Figure 5 shows two remote sensing images taken before and after the 2010 Haiti earthquake. (a) is the initial matching result, and (b) is the matching result after applying the TTM algorithm.

Applying the above algorithm to remote sensing images, we divide the experimental images into two groups: 1) remote sensing images of different spectra; 2) remote sensing images before and after disasters. Finally, the Precision, recall and RMSE indicators are used to evaluate the results, and compared with the classic RANSAC, GTM and Bi-SOGC methods.

Example 1: Different spectral remote sensing images

15 pairs of remote sensing images with different spectra are used, and the error matching ratio is set at 5%-95%. Matching examples and accuracy results are shown in Figure 7 and Figure 8. In Figure 7(b), (c), (d), and (e), the solid line connection is the correct match, and the dotted line connection is the wrong match. By comparing the present invention with other methods, it can be seen that the present method has achieved the highest precision in terms of Precision, recall and RMSE.

Example 2: Remote sensing images before and after disasters

15 pairs of remote sensing images before and after the disaster are used, and the error matching ratio is set at 5%-95%. Matching examples and accuracy results are shown in Figure 9 and Figure 10. In Figure 9(b)(c)(d)(e), the solid line connection is correct matching, and the dotted line connection is wrong matching.

Through comparative experiments, it can be proved that the method of the present invention has the best performance in terms of accuracy index, especially when the proportion of correct matching is small.

Claims (2)

1. a point matching method of remote sensing image, is characterized in that: the method comprises the steps: Step (1): The initial matching point set is two corresponding point sets extracted from the input in the reference image and the sensing image, in the initial matching point set P={p _i } and P'={p _i '} , i=1,2,3...N, N is the number of points in the point set, three point pairs are selected arbitrarily to form many corresponding triangle pairs, and the three point pairs selected arbitrarily are expressed as (p _i, p _j, p _k ) and (p _i ', p _j ', p _k '); Step (2): Calculate the corresponding angle difference and the corresponding side length ratio difference between the two triangles formed: <mrow><mfencedopen='{'close=''><mtable><mtr><mtd><mi>&amp;Delta;</mi><msub><mi>&amp;theta;</mi><mn>1</mn></msub><mo>=</mo><mo>|</mo><mo>&amp;angle;</mo><msub><mi>p</mi><mi>j</mi></msub><msub><mi>p</mi><mi>i</mi></msub><msub><mi>p</mi><mi>k</mi></msub><mo>-</mo><mo>&amp;angle;</mo><msup><msub><mi>p</mi><mi>j</mi></mi>msub><mo>'</mo></msup><msup><msub><mi>p</mi><mi>i</mi></msub><mo>'</mo></msup><msup><msub><mi>p</mi><mi>k</mi></msub><mo>'</mo></msup><mo>|</mo></mtd></mtr><mtr><mtd><mi>&amp;Delta;</mi><msub><mi>&amp;theta;</mi><mn>2</mn></msub><mo>=</mo><mo>|</mo><mo>&amp;angle;</mo><msub><mi>p</mi><mi>i</mi></msub><msub><mi>p</mi><mi>j</mi></msub><msub><mi>p</mi><mi>k</mi></msub><mo>-</mo><mo>&amp;angle;</mo><msup><msub><mi>p</mi><mi>i</mi></msub><mo>'</mo></msup><msup><msub><mi>p</mi><mi>j</mi></msub><mo>'</mo></msup><msup><msub><mi>p</mi><mi>k</mi></msub><mo>'</mo></msup><mo>|</mo></mtd></mtr><mtr><mtd><mi>&amp;Delta;</mi><msub><mi>&amp;theta;</mi><mn>3</mn></msub><mo>=</mo><mo>|</mo><mo>&amp;angle;</mo><msub><mi>p</mi><mi>i</mi></msub><msub><mi>p</mi><mi>k</mi></msub><msub><mi>p</mi><mi>j</mi></msub><mo>-</mo><mo>&amp;angle;</mo><msup><msub><mi>p</mi><mi>i</mi></msub><mo>'</mo></msup><msup><msub><mi>p</mi><mi>k</mi></msub><mo>'</mo></msup><msup><msub><mi>p</mi><mi>j</mi></msub><mo>'</mo></msup><mo>|</mo></mtd></mtr></mtable></mfenced><mo>-</mo><mo>-</mo><mo>-</mo><mrow><mo>(</mo><mn>1</mn><mo>)</mo></mrow></mrow> <mrow><mfencedopen='{'close=''><mtable><mtr><mtd><mi>&amp;Delta;</mi><msub><mi>l</mi><mn>1</mn></msub><mo>=</mo><mo>|</mo><msub><mi>p</mi><mi>i</mi></msub><msub><mi>p</mi><mi>k</mi></msub><mo>/</mo><msup><msub><mi>p</mi><mi>i</mi></msub><mo>'</mo></msup><msup><msub><mi>p</mi><mi>k</mi></msub><mo>'</mo></msup><mo>-</mo><msub><mi>p</mi><mi>i</mi></msub><msub><mi>p</mi><mi>j</mi></msub><mo>/</mo><msup><msub><mi>p</mi><mi>i</mi></msub><mo>'</mo></msup><msup><msub><mi>p</mi><mi>j</mi></msub><mo>'</mo></msup><mo>|</mo></mtd></mtr><mtr><mtd><mi>&amp;Delta;</mi><msub><mi>l</mi><mn>2</mn></msub><mo>=</mo><mo>|</mo><msub><mi>p</mi><mi>i</mi></msub><msub><mi>p</mi><mi>k</mi></msub><mo>/</mo><msup><msub><mi>p</mi><mi>i</mi></msub><mo>'</mo></msup><msup><msub><mi>p</mi><mi>k</mi></msub><mo>'</mo></msup><mo>-</mo><msub><mi>p</mi><mi>k</mi></msub><msub><mi>p</mi><mi>j</mi></msub><mo>/</mo><msup><msub><mi>p</mi><mi>k</mi></msub><mo>'</mo></msup><msup><msub><mi>p</mi><mi>j</mi></msub><mo>'</mo></msup><mo>|</mo></mtd></mtr><mtr><mtd><mi>&amp;Delta;</mi><msub><mi>l</mi><mn>3</mn></msub><mo>=</mo><mo>|</mo><msub><mi>p</mi><mi>j</mi></msub><msub><mi>p</mi><mi>k</mi></msub><mo>/</mo><msup><msub><mi>p</mi><mi>j</mi></msub><mo>'</mo></msup><msup><msub><mi>p</mi><mi>k</mi></msub><mo>'</mo></msup><mo>-</mo><msub><mi>p</mi><mi>i</mi></msub><msub><mi>p</mi><mi>j</mi></msub><mo>/</mo><msup><msub><mi>p</mi><mi>i</mi></msub><mo>'</mo></msup><msup><msub><mi>p</mi><mi>j</mi></msub><mo>'</mo></msup><mo>|</mo></mtd></mtr></mtable></mfenced><mo>-</mo><mo>-</mo><mo>-</mo><mrow><mo>(</mo><mn>2</mn><mo>)</mo></mrow></mrow> Step (3): and Represent the horizontal and vertical coordinates of points p _i and p _i ', and set the reference in the two triangles respectively Point p _d and p _d ', the coordinate calculation formula is as follows: <mrow><mfencedopen='{'close=''><mtable><mtr><mtd><msub><mi>x</mi><msub><mi>p</mi><mi>d</mi></msub></msub><mo>=</mo><mfrac><mn>1</mn><mn>2</mn></mfrac><mrow><mo>(</mo><msub><mi>x</mi><msub><mi>p</mi><mi>i</mi></msub></msub><mo>+</mo><msub><mi>x</mi><msub><mi>p</mi><mi>j</mi></msub></msub><mo>)</mo></mrow></mtd></mtr><mtr><mtd><msub><mi>y</mi><msub><mi>p</mi><mi>d</mi></msub></msub>msub><mo>=</mo><mfrac><mn>1</mn><mn>2</mn></mfrac><mrow><mo>(</mo><msub><mi>y</mi><msub><mi>p</mi><mi>i</mi></msub></msub><mo>+</mo><msub><mi>y</mi><msub><mi>p</mi><mi>j</mi></msub></msub><mo>)</mo></mrow></mtd></mtr></mtable></mfenced><mo>-</mo><mo>-</mo><mo>-</mo><mrow><mo>(</mo><mn>3</mn><mo>)</mo></mrow></mrow> <mrow><mfencedopen='{'close=''><mtable><mtr><mtd><msub><mi>x</mi><msup><msub><mi>p</mi><mi>d</mi></msub><mo>'</mo></msup></msub><mo>=</mo><mfrac><mn>1</mn><mn>2</mn></mfrac><mrow><mo>(</mo><msub><mi>x</mi><msup><msub><mi>p</mi><mi>i</mi></msub><mo>'</mo></msup></msub><mo>+</mo><msub><mi>x</mi><msup><msub><mi>p</mi><mi>j</mi></msub><mo>'</mo></msup></msub><mo>)</mo></mrow></mtd></mtr><mtr><mtd><msub><mi>y</mi><msup><msub><mi>p</mi><mi>d</mi></msub><mo>'</mo></msup></msub><mo>=</mo><mfrac><mn>1</mn><mn>2</mn></mfrac><mrow><mo>(</mo><msub><mi>y</mi><msup><msub><mi>p</mi><mi>i</mi></msub><mo>'</mo></msup></msub><mo>+</mo><msub><mi>y</mi><msup><msub><mi>p</mi><mi>j</mi></msub><mo>'</mo></msup></msub><mo>)</mo></mrow></mtd></mtr></mtable></mfenced><mo>-</mo><mo>-</mo><mo>-</mo><mrow><mo>(</mo><mn>4</mn><mo>)</mo></mrow></mrow> Step (4): Set triangle similarity threshold θ _b = 5°, l _b = 0.1, and then according to the order consistency and conditions of the three vertices relative to the reference point △θ ₁ , △θ ₂ , △θ ₃ ∈[0 ,θ _b ], △l ₁ ,△l ₂ ,△l ₃ ∈[0,l _b ] to filter out the qualified triangle pairs; Step (5): In the selected triangle pair, count the number of occurrences of each point pair and establish a frequency histogram; Step (6): Set the threshold k _p according to the experimental analysis, and remove the point pairs whose occurrence times are less than k _p ; Step (7): Calculate the model conversion parameters according to all point pairs at this time, the point pairs are expressed as P _o and P _o ', and the number is n; T=argminE(T(θ)) (10) in, Represents the root mean square error, T(θ) represents the conversion parameters obtained by the point set P _o and P _o 'using the least square method, d(i)=||T(θ)(p _i ')-p _i ||, i=1,2,3,...n represents the conversion error of each point pair; Step (8): Set the threshold E _t , when E(T(θ))>E _t , the point pair with the largest d(i) value will be removed, and the point set will be updated at the same time, when the condition E(T (θ))≤E _t , the iteration ends, and two correct matching point sets are obtained at this time. 2. The point matching method of a remote sensing image according to claim 1, characterized in that: the method is used to screen out the correct point matching set from the initial matching point set of the remote sensing image with a large background change.