Subpixel matching method for remote sensing image of ground features based on geographic information

2.1.1 Low-frequency suppression

Because most of the ability of remote sensing image is concentrated in the low-frequency part (the center of the spectrum image), the amplitude at the low-frequency part is larger than that at other spectrum positions in the spectrum image of remote sensing image. Let f be the original remote sensing image, F be the spectrum image obtained by Fourier transform, and |F| be the amplitude of F. Assuming that the point (x, y) is the coordinate of the maximum value in |F|, that is, the coordinate at the low frequency, and the suppression function is taken as follows:

where 0 < k < 1 is the suppression constant. A smaller k indicates a greater degree of low-frequency suppression and vice versa. If k is too large (close to 1), the enhancement effect is not obvious; if k is too small (close to 0), the enhancement result will be distorted. The value of k can be determined according to the modulation transfer function (MTF) value at the Nyquist frequency of the imaging system. When the MTF value at the Nyquist frequency is small, it means that the loss of high-frequency capability is serious after image degradation, then k should take the smaller value, and vice versa [7].

2.1.2 Gamma correction

The frequency spectrum image F′ after low-frequency suppression is transformed by inverse Fourier transform to obtain the low-frequency suppression image g, and then the final enhanced image r is obtained by Gamma correction of g as follows:

where g′ is the normalized low-frequency suppressed image and γ is the Gamma correction factor. To keep the energy of the enhanced image and the energy of the original image suppressed, take

In the formula, m _g′ and m _f are gray mean values of g′ and f, respectively, and τ is constant. The constant τ is used to adjust the overall brightness of the image. If the overall brightness of the original remote sensing image is dimmed due to the decrease of the remote sensor gain, τ is taken as a negative constant to make γ close to 0, to improve the overall brightness of the enhanced remote sensing image r. Assuming the value of τ is 0, the Gamma correction factor γ > 1, the mean gray value of the normalized low-frequency suppression image is greater than the mean gray value of the original remote sensing image. After the correction of Gamma correction factor, the overall gray value of g′ will be reduced, approaching the mean gray value of the original remote sensing image, to keep the energy unchanged before and after the enhancement as much as possible; when the mean gray value of the normalized low-frequency suppression image is small than that of the original remote sensing image, γ < 1. After Gamma correction, the overall gray value of g′ will be improved, approaching the gray value of the original remote sensing image [8].

2.2.1 Boundary geographic information extraction in multisensor image

In this paper, an improved thresholding algorithm based on histogram exponential convex hull is used to automatically extract the boundary geographic information in the enhanced remote sensing image r. The boundary geographic information extracted in this paper is mainly the boundary information of coastline, lake, and other remote sensing images with obvious geographical indications [9].

Let h(k), k ∈ [0, n] be the histogram function of the enhanced remote sensing image r, n be the maximum gray level, and h ¯ ( k ) be the convex hull of the polygon (that is, the minimum convex polygon containing h(k)). On the basis of polygon convex hull h ¯ ( k ) , the definition of histogram exponential convex hull h ¯ e ( k ) can be derived. The process of getting histogram exponential convex hull h ¯ e ( k ) is as follows, that is to get the polygon convex hull In h ¯ ( k ) of In h(k) first, and then get h ¯ e ( k ) by exp [ In h ¯ ( k ) ] . Assuming that the minimum and maximum gray levels of h(k) ≠ 0 are L _min and L _max, respectively, for h(k), k ∈ [L _i , L _i+1] and L _min ≤ L _i ≤ L _i+1 ≤ L _max, the index convex h ¯ e ( k ; L i , L i + 1 ) within the gray level range [L _i , L _i+1] can be calculated as follows:

Let the residual of the exponential convex hull in the gray level range [L _i , L _i+1] be

It can be used to describe the approximation degree of the exponential convex hull h ¯ e ( k ; L i , L i + 1 ) in the gray level range [L _i , L _i+1] to the histogram h(k) of the enhanced remote sensing image r, and it also reflects the convex and concave degrees of h(k) in the gray level range [L _i , L _i+1]. Compared with the polygon convex hull h ¯ e ( k ; L i , L i + 1 ) , the exponential convex hull h ¯ e ( k ; L i , L i + 1 ) can better approach the original histogram h(k), which is conducive to the threshold selection of thresholding segmentation. Therefore, this paper chooses the thresholding segmentation algorithm based on the convex hull selection threshold of histogram exponential.

In the thresholding segmentation algorithm based on histogram exponential convex hull, it is assumed that the gray level range [L _min, L _max] is divided into N categories, the threshold thresholds selected are t ₁, t ₂,...,t _N−1, respectively, and the overall residual of the exponential convex hull is defined as follows:

Then, the threshold t ₁, t ₂,...,t _N−1 selected shall be such that

In the boundary geographic information extraction of this paper, the enhanced remote sensing image r needs to be divided into two categories (i.e., N = 2), so as long as the threshold t ₁ is determined, the following can be achieved:

Then, the image can be divided into two categories by threshold t ₁, based on which boundary geographic information can be further extracted.

In the process of thresholding the actual t ₁ image, Landsat MSS and TM image, we find that the above method of thresholding by histogram exponential convex hull is often very sensitive to noise; sometimes, it will lead to improper selection of threshold, making the image on histogram almost fall on one side of thresholding threshold, while the other side of the threshold is almost none. Therefore, histogram must be preprocessed before threshold segmentation [2].

Therefore, the concept of population is introduced, that is to eliminate the influence by population test. The population P _i at the gray value i is defined as

Population P _i represents the number of pixels on both sides of a certain threshold in the remote sensing image r. Before getting the histogram exponential convex hull, a part of gray level pixels with a very small population are removed, and then the histogram is balanced accordingly. The advantage of this method is that it can eliminate the influence of noise in the high brightness area (larger gray value) and low brightness area (smaller gray value) of the image, so that the selected threshold does not fall in the false area where the noise is located.

2.2.2 Boundary image matching

Boundary image matching is a kind of geographical boundary matching in remote sensing images. The basic principle of remote sensing image matching is to find out the location with the largest similarity measurement value or the smallest difference between the two remote sensing images (measured map, benchmark map) according to the similarity measurement of geographical information in the two remote sensing images, while the common similarity measurement includes the minimum distance measurement, the maximum correlation measure, and invariant moment measure, etc. According to the characteristics that the boundary image obtained after threshold segmentation is a binary image, this paper adopts a simple and intuitive similarity measurement-pair function measurement method, which has been applied to the automatic matching of boundary image and achieved good results [10]. Compared with several common similarity measurement methods, the function measurement method is more suitable for binary boundary image matching and has the advantages of simple, fast, and strong denoising performance. This method is introduced as follows:

If the measured image size of the remote sensing image is N ₁ × N ₂, the reference image size is M ₁ × M ₂ (usually M ₁ > N ₁, M ₂ > N ₂), and the gray values at the coordinates (x, y) of the measured image and the reference image are g(x, y) and f(x, y) respectively, then the measured image can be expressed as the N ₁ × N ₂ × 1-dimensional image vector, recorded as G, i.e. G = [g(0, 0), g(0, 1),..., g(1, 0), g(1, 1),..., g(N ₁ − 1, N ₂ − 1)]; similarly, the starting position is set as (u, v) and the reference subimage with the same size as the measured image is also expressed as N ₁ × N ₂ × 1, i.e., F _u,v .

If the gray level range of the image element in the measured remote sensing image and the reference subimage is [0, n], then the corresponding pair function N _ws (u, v) can be established between any gray level w in the measured image and any gray level s in the reference subimage whose starting position is (u, v). The pair function N _ws (u, v) is defined as the number of image element composed of pixels with gray level w in the measured image and pixels with gray level s at corresponding points in the reference subimage with starting position (u, v), that is,

where N w s ( ∘ ) is the operation to calculate the number of set elements in brackets, w, s ∈ [0,n].

If the measured image and the reference image are binary boundary images. Then, the definition of the function can refer to the definition of gray image, which is recorded as N _ws (u, v), w, s ∈ [0,1], where w, s = 1 represents boundary pixel.

Several kinds of measurement algorithms can be defined by the pair function. In this paper, the product measurement algorithm of the pair function is adopted. For a binary image, the product measurement algorithm of the pair function can be expressed as follows:

where R(u, v) is the product of the ratio of the number of matched pixels’ pairs to the number of pairs of matched pixels. The larger R(u, v) is, the better the similarity is.

When the product measure algorithm of function is used for binary boundary image matching, if there are two or more test points with the same maximum correlation coefficient R(u, v) in the reference image, this paper adopts the method of measuring the boundary center to find the best matching position (u′, v′), that is, the obtained measured image and the boundary barycenter position ( x ¯ g , y ¯ g ) and ( x ¯ f , y ¯ f ) of the reference subimage in the test point (u, v) are as follows:

Then, the best matching point (u′, v′) is obtained by the following method, that is, the matching position D(u, v), min u , v ( u , v ) :

2.2.3 Image matching and positioning based on least squares

After the boundary image matching based on geographic information, we can get a series of image block’s center coordinates (x _i , y _i ) and (X _i , Z _i ) corresponding to the measured image and the reference subimage that achieve the best matching under the function product measure. They will be used as the control block pairs between the measured image and the reference image. According to the corresponding geometric transformation relationship between the measured image and the reference image established above, the matching and positioning of the measured remote sensing image and the reference image can be carried out under the least square error criterion, also known as the least square matching and positioning [11,12]. It is assumed that the corresponding geometric transformation relationship between the measured image and the reference image is expressed as a bivariate quadratic polynomial, namely,

where a and b represent subpixels. The problem of matching and positioning between the measured image and the reference image can be transformed into the sum of the linear least-squares solutions A and B of the equations { W A = U W B = V , where

Then, there are the least-squares solutions of A and B

where U and V represent pixels. It should be pointed out that in the least-square image matching and positioning in this paper, it is assumed that the corresponding geometric transformation relationship between the measured image and the reference image can be approximated by bivariate quadratic and bivariate first-order polynomials [13,14,15,16]. The basic starting point is that the optional correction polynomials (bivariate quadratic and bivariate first-order polynomials) are used to directly simulate the geometric distortion model of the measured image. It is considered that the geometric distortion of the measured remote sensing image can be regarded as the comprehensive effect of translation, scaling, rotation, affine, skew, bending, and the basic deformation of a quadratic polynomial [17,18,19,20,21,22]. In this way, we can avoid the complex imaging spatial geometric process of the measured remote sensing image to establish the corresponding geometric transformation relationship between the measured image and the reference image, to reduce the mismatch and improve the subpixel matching accuracy of the image [23,24,25,26].