COMBINED MULTI-VIEW MATCHING ALGORITHM WITH LONG-STRIPS OF SATELLITE IMAGERY FROM DIFFERENT ORBITS

Existing matching algorithms aim to match conjugate points among overlapping satellite scenes acquired from the same orbit and can generally achieve good matching performance. Unfortunately, no algorithm can avoid the difficulty of simultaneously processing the data sets of long-strip imagery acquired from different orbits. In this paper, the combined matching algorithm we propose introduces the LBP/C operator, which, when combined with existing feature detectors for the first time, can make possible the extraction of more stable interest points and candidates. At the same time, based on the typical characteristics of Chinese satellite imagery, we improved the filter method and achieved an effective combination of several image matching algorithms. A comparison among several kinds of matching transfer modes was presented; and to evaluate this algorithm, Chinese Mapping Satellite-I data are used as the reference data.


INTRODUCTION
Image matching was a key technology related to many problems in computer vision and photogrammetry, including object recognition, motion tracking, digital surface model (DSM) generation, etc. (Lowe, 2004).In the past 20 years, research related to processing satellite imagery has become a key focus.Due to the position and attitude information under fixed sampling intervals, the current matching algorithms are focused on area-based algorithms, especially the approximate epipolar geometric constraints (Heipke, 1996).These existing algorithms can achieve satisfactory matching results in the case of overlapping scenes acquired from the same orbit, but they cannot solve the following problems: 1) how to match with multiple long-strip satellite imagery simultaneously 2) how to eliminate the radiometric and geometric differences between adjacent orbits and achieve high-accuracy results.Moreover, the typical characteristics of imagery acquired from Chinese satellites (e.g., ambiguity in smooth areas, low direct georeferencing accuracy, and obvious differences in contrast) bring more challenges to image matching.Consequently, finding solutions for the above-mentioned problems in photogrammetry was an urgent issue.In this study, we developed a practical matching algorithm based on the approaches in the literature.First, we divided the reference strip of imagery in each orbit into independently overlapping image patches; then, each patch was filtered to enhance the existing texture patterns and improve the noisesignal ratio (Baltsavias, 1991).After the generation of image pyramids, we combined the Harris-Z detector (Bellavia, 2008) with the LBP/C operator (Timo, 2002) to extract stable interest points.At the same time, the Förstner interest operator (Förstner, 1986) was used to refine these interest points.Then, through the global digital elevation model (Global DEM), the approximate geography scope covered by each image patch was computed.Moreover, the Geometrically Constrained Cross-Correlation (GC 3 ) algorithm (Zhang, 2005 and2006) was extended to match multiple long-strip imagery; and finally, several strategies were proposed to detect and eliminate the mismatches.By comparing the results of several matching transfer modes, we determined the optimal manner in which to optimize the matching speed and reduce the number of mismatches.
In the second part, the new matching algorithm is described in detail, and the validity of the proposed matching algorithm is verified with Chinese Mapping Satellite-I data.

METHODOLOGY
In this paper, the proposed algorithm includes four parts: 1) Image Pre-Processing; 2) Image Matching Procedure; 3) Mismatch Detection and Elimination.

Image Pre-Processing
Image pre-processing method includes two parts: 1) Modified Wallis Filter; 2) Interest Points and Candidates Detection and Extraction

Modified Wallis Filter
Compared with aerial imagery and close-range imagery, the characteristics of satellite imagery generally include the following: 1) low spatial resolution, 2) unsymmetrical contrast of light and shade, 3) fuzzy texture structure.These factors seriously affect the reliability of feature extraction and image matching.
To solve these problems mentioned above, the Wallis filter is introduced.It can force the mean and contrast of an image to fit some given values effectively.Baltsavias (1991) used this filter to enhance the image texture pattern and eliminate radiometric difference.The general form of the Wallis filter is given by: Where In Equation 1, x y w ( , ) However, the existing Wallis filter has some problems: 1) the edges of roofs tend to expand, which causes deviations from their geometrical position; and 2) the image noise is partially enhanced, which hurts the positioning accuracy of the features extracted and has a detrimental effect on image matching.
We found that the additive parameter 0 r causes the above problems.Owing to its restrictions on the gray level, unreasonable selection of 0 r causes the excessive loss of gray information for the following reasons.First, 0 r is mainly dependent on parameter b ; secondly, the multiplicative parameter 1 r can change image contrast, whose value is dependent on parameter c ; and lastly, the high texture fidelity and enhance image contrast must be kept under the filter's capacity.In this study, we used the following method to improve the Wallis filter.
A highly adaptive parameter determination strategy: In the case of an image f(x,y), we classified all pixels of this image into two sets (C1,C2).2002) developed the LBP operator, which allows for detecting "uniform" patterns in circular neighborhoods of any quantization of the angular space and spatial resolution.Moreover, he enhanced the performance of the LBP operator through combining it with a rotation invariant local contrast (LC) measure that characterizes the contrast of local image texture.
In the case of an image, texture T was the joint distribution of the gray levels: In fact, without losing information, the signed differences between the gray value of c g and P g (P=0,…P-1) were not affected by variation in the mean luminance.Hence, the joint difference distribution was invariant against the grayscale shifts.By considering the signs of the differences, instead of their exact values, the operator could achieve the scale invariance: Where By assigning a binomial factor 2 P for each sign s , we transformed Equation 9 into a unique LBP value.In this study, we set (P=12,R=1.5): 11) Ojala designated patterns that have U values of at most 2 as "uniform" and proposed the following operator for grayscale and rotation invariant texture description instead of The 2 , riu P R

LBP
operator was a grayscale invariant measure.It was an excellent measure of the spatial pattern, but it discarded contrast.So a rotation invariant measure was introduced to characterize the contrast of the local image texture: , was expected to be a very powerful rotation invariant measure of local image texture.When we compared feature detectors, Harris-Z was found to be more suitable to detect the interest points for large data sets of satellite imagery in our study.We also introduced the Förstner detector at this point and the surface fitting approach to get sub-pixel positioning.In this study, we combined Harris-Z with LBP/C to detect the stable interest points and candidates that optimize the validity of matching, which included the following five steps: 1) To ensure uniform distribution of the interest points, we used the grid of the feature extraction as a region unit, and a fixed number of features were extracted in each unit independently.
2) The interest strengths and LBP/C values in a region unit were calculated by the Harris-Z detector and the LBP/C operator.
3) Some points were saved when their LBP/C values were in a certain range.Then, we chose the point with the local maximum interest strength as the interest point that corresponded to the region unit.As a candidate, we chose the point with the maximum interest strength without consideration for its LBP/C value.4) We used the Förstner detector and the surface fitting approach to obtain the sub-pixel positioning of the interest point and candidate.5) Once the interest point failed to match or the accuracy of this match was low, the candidate was matched rather than the interest point; and the match with the higher accuracy was saved.

Image Matching Procedure
Due to linear pushbroom imaging mode, we can calculate the exterior orientation corresponding to every scanline through lagrangian interpolation.Thereby, in this paper, the proposed matching algorithm is based on the geometrically constrained cross-correlation.The algorithm includes the following several studies:

Automatic Computation for Approximate Elevation Range and Elevation
Step through Global DEM According to the position and attitude information of the scan lines, we use a hypothesis height to the frontal-project four corners of the corresponding reference image patch onto the ground.The approximate geography scope covered by the image patch was determined, and the maximum elevation H max and the minimum elevation H min are achieved through retrieval from Global DEM.
Then, we frontal-projected an interest point extracted onto the ground using H max and H min , and back-projected these object space points onto all of the search strips.
It was an iterative process: once i H pitch was an invalid value, we had to repeat the above procedure using another interest point.

Geometric and Radiometric Distortion Procedure
There were geometric and radiometric distortions in the case of strips acquired from different camera lenses.When the correlation window is defined in the reference image, this kind of distortion leads to an irregular and discontinuous shape of the corresponding window in the search images, which cannot be matched in a straightforward manner.Gruen and Baltsavias (1985) extended the LSM method, which compensates this distortion and relaxes the strict geometric and radiometric assumption for normal cross-correlation methods.So, we introduced this method to eliminate the window warping: 1) In the reference image, we defined a correlated window Л and a search window Г in which an interest point or candidate was central, then the corresponding pixel coordinates of the four corners and the center of Г were achieved by projection in the search image, which formed an irregular polygon Г', as shown in Fig. 2. Accordingly, we calculated the parameters of the affine geometric transformation and linear radiometric transformation between these images.
Where 0 a , 1 a , 2 a , 0 b , 1 b , 2 b were the parameters of the affine geometric transformation; 0 h , 1 h were the parameters of the linear radiometric transformation.
2) Through Equations 17 and 18, the relationship of each pixel between Г and Г' was determined.The gray values in Г' were interpolated by bilinear interpolation, and these gray values, recalculated by Equation 19, were assigned to the corresponding pixel in Г.  3) Once the angle θ was greater than 90 degrees on the very steep terrain surface, as shown in Fig. 3, W' became zero and mafe the procedure impossible.In this case, we found that when we iteratively computed the parameters of the transformation, the convergent speed was slow, and even the computation was divergent.So once the threshold of iteration times was crossed, we considered θ≥90 and excluded it.4) The correlation window and search window are both square as shown in Fig. 4. Therefore, we could do a straightforward match and reduce the computational cost.

Extended Geometrically Constrained Cross-Correlation Strategy
Based on the area-based matching techniques, Zhang developed GC 3 (Geometrically Con-strained Cross-Correlation) algorithm.With the GC 3 matching procedure, multiple images can be matched simultaneously and the epipolar constraint is integrated implicitly.As a similarity measure, SNCC (Sum of Normalized Cross-Correlation) greatly decreases mismatches caused by approximate texture distinction and occlusion problem.However, this algorithm causes high computational costs in determining the correct height value and has redundance for its quality measuring procedure.Basically we extended this algorithm and save the computational costs: 1) For each reference image patch, the approximate elevation range and elevation step can be determined.After preprocessing the reference image patch, for each interest point or candidate, the corresponding approximate epipolar lines in the search strips can be further determined.
2) According to the position of approxi-mate epipolar lines, we define the search image patches, and pre-process them.
3) Based on geometric and radiometric distortion procedure, we calculate the parameters of affine geometric transformation and linear radiometric transformation between the defined windows.4) We combine a coarse-to-fine hierarchical approach with the efficient implementation of the epipolar geometry to match each point on the approximate epipolar lines.Once normalized cross-correlation coefficient in the certain level of pyramids is lower than threshold, we repeat the above procedures with different elevation in the certain range.5) Until accomplishing the matching procedure, we compute the SNCC of all matching candidates.6) Determine the sites of local maxima of the SNCC function and fit a smooth quadratic function for each local maximum with a local neighborhood of that maximum value.7) If the second SNCC peak is less than half or one third of that of the first SNCC peak, the peak that have the largest SNCC value should represent the correct match.If not, we match inversely and if the difference between the two matches is less than 1.5 pixels in image space, the candidate should be the correct match.8) Repeat 1)-7), until all interest points and candidate are matched.

Mismatches Detection and Elimination
After finishing the image matching procedure, there are inevitably some blunders that do not conform to the geometrical relationship in an object space, which we call mismatches.To maintain the matching accuracy, we must eliminate mismatches as accurately as possible.Thus, we classify all of the matches into two types, interior matches and adjacent matches.In this section, we present two strategies to eliminate mismatches

Mismatches Elimination Strategy within Orbit:
In this study, we introduce the multi-rays forward intersection principle to eliminate mismatches and realize the object fusion for the matching results.The implementation strategy was as follows: 1) If an interest point or candidate finds multiple correspondences on the search strips, we assume that the exterior orientation elements corresponding to each image are known.We take the 3D coordinates of all matches in the object space as unknown; and based on the collinearity equations, we then list the error equations: , y l are differences between observations and approximations.
2) If the number of the correspondences is n, we can list (2n+2) error equations.To solve the iterative computation for unknown numbers, the least square adjustment was introduced.In the iterative procedure, we set the weights of each observation of the image points through a hypothesis test of the posterior variance component.If there are grosses in the candidates, the observation weights of the candidates become smaller and smaller until 0 was reached, which makes no sense in the adjustment process in order to eliminate the effect of a false match on the calculation of effectiveness.3) After iterative convergence was achieved or the iterative times exceed the limited value, the root mean square (RMS) value was calculated.Once the RMS exceeds the given threshold, the corresponding matches are eliminated.

Mismatches Elimination Strategy between Adjacent Orbits
Once the position and attitude parameters acquired from different orbits exhibited obvious system shift, the matches between the adjacent orbits could not be intersected together in the object space.If we used the above strategy, we would face the problem of iterative computation convergence.Therefore, we introduced an Object Space Shift (OSS) strategy to eliminate mismatches between the adjacent orbits: 1) The correspondences P1, P2, …, Pn between the adjacent orbits were intersected within each orbit, and the shift of each point in the object space were calculated ( X were calculated in Equation 22: The mean square values in each direction x m , y m , z m were calculated in Equation 23, and the RMS error m between the adjacent orbits were calculated in Equation 24. 2) Once the shift i T was more than triple error m , we deemed this correspondence as a mismatch.

Test Data Description
Figure 4 shows the real data for experiment, which consists of two long strips of imagery acquired by Chinese Mapping Satellite.These data are imaged on October, 2010, which covers Harbin in China, the coverage area is about 25000km 2 .The satellite sensor works in pushbroom mode, and the stereo module consists of three lenses with one CCD line scanner for each, which provide a forward-, a nadir-and a backwardlooking view.The degree of overlapping within orbit is about 98%, and the degree between adjacent orbits is about 10%.The parameters of CCD line sensors are shown in table 1.For testing and verifying the accuracy and reliability of the proposed algorithm, the geographical region covered by the test data includes not only farms and forests which have strong texture repeatability, but also residential quarters, mountain countries and smoothly open areas which are known as difficult terrain for matching.In addition, in order to analyze the matching accuracy quantitatively, the test data includes 33 planar control points and 15 elevation control points.After finishing the matching procedure, the ground control points are measured interactively on the test imagery, and then the bundle block adjustment is performed with matching points and control points.

Image Matching Testing
In order to improve the performance of the proposed matching algorithm, we compare different parameters in two ways: the number of pyramid levels and the size of the relative window at the top level.Figure 5 shows the results in which the number of pyramid levels and the size of the relative window at the top level are varied.In graph (a), when the number is 3, the successful rate of matching is the highest toward the different terrain types.In graph (b), it is poor at the successful rate when the relative windows is set 7×7 pixels, but the results are not improved till 13×13 pixels.After that, a larger size of relative window can actually hurt the matching.Although the successful rate for 13×13 is a little lower than 11×11 and 9×9 on mountain countries and farms respectively, the overall successful rate is optimum.Therefore, through this experiment, we generate three pyramid levels and set the size of the relative window R×R as 13×13.Moreover, the relative window is shrank to (R-2)×(R-2) level to level, and the corresponding search window is (2×R-1)*(2×R-1).The matching results are listed in table 3. Table 3 illustrates the statistics of image matching for the test date on five types of terrain, through analysis for the matching results, we can see that: 1) For farms and forests which have strong texture repeatability, the range of the successful rate of matching is from 76.9% to 90.7%.Especially in areas covered by farms, the RMS error is within 0.5 pixel, which means that the proposed algorithm can avoid mismatching problem caused by texture repeatability.
2) For residential quarter and smoothly open areas, there are a few false matches, but most of matches are fine, the range of the successful rate is from 77.6% to 94.7%, the range of RMS error is from 0.21 to 0.47 pixel.Thus, the matching result is satisfactory.
3) For mountain countries, because of the rise or descend on terrain greatly, there are some obviously false matches, which greatly affect the precision of matching.But the RMS error is still 0.74 pixel, which proves that the matching algorithm can be applied to this terrain.

CONCLUSION
In this study, we develop a practical algorithm for matching with large data sets of long-strips of satellite imagery.Through some experiments using real data set from Chinese Mapping Satellite, the proposed algorithm are proved to be valid, and the satisfactory accuracy is achieved under fully automatic computation.
At the present stage, the proposed methods mostly focus on multi-view matching with long strips of satellite imagery acquired from different orbits.In the futhre, we will keep on this research and realize the matching algorithm that support various satellite sensors, even various imaging modes.
the mean and standard deviation, c was the contrast expansion constant, and b was the brightness forcing constant.

c
means the set of high-pass.Then, b and c were calculated by the following equationtool for texture classification, the Local Binary Pattern (LBP) operator was an excellent measure of the spatial structure of local image texture.Ojala (

Figure 2 .
Figure 2. The principle of the affine geometric transformation and linear radiometric transformation between the reference image and search image

Figure 3 .
Figure 3. Image selection or exclusion based on the terrain slope and the sensor geometry.The size of the corresponding correlation window in the search image depends on the angle θ between the normal vector of the local object surface and the corresponding imaging ray.
related to the line elements of the exterior orientation; x l

Figure 4 .
Figure 4.The Nadir-looking strips in the adjacent orbits: (a) The nadir-looking strip acquired from the first orbit; (b) The nadir-looking strip acquired from the second orbit; (c) The enlarged drawing of the overlapping area between the two strips

Figure 5 .
Figure 5.The successful rate of matching using the number of pyramid levels and the size of the relative window at the top level: (a) the effect of different number of pyramid levels; (b) the effect of different size of relative window at top level

Table 1 .
Camera parametres of the three linear array sensor for Chinese Mapping Satellite

Table 3 .
Matching accuracy quantitative evaluation for different types of terrain