JP2003141507A - Precise geometric correction method for landsat tm image and precise geometric correction method for satellite image - Google Patents

Abstract

PROBLEM TO BE SOLVED: To provide a precise geometric correction method for Landsat TM images enabling accurate geometric correction for an entire area of interest. SOLUTION: A geometric transformation containing a translation parameter included in system information provided along with a Landsat TM image is as follows: p-p0 =[cosθ0 *(x-x0 )-sinθ0 *(y-y0 )]/Δ; l-l0 =[-sinθ0 *(x-x0 )-cosθ0 *(y-y0 )]/Δ. A solar direct illumination image created using a numerical altitude model is used as a simulation image. An appropriate translation parameter is determined while using as an evaluation function a correlation function of the pixel values of the simulation image and the Landsat TM image.

Description

Detailed Description of the Invention

[0001]

TECHNICAL FIELD The present invention relates to Landsat TM.
The present invention relates to an image precision geometric correction method and a satellite image precision geometric correction method.

[0002]

2. Description of the Related Art Landsat TM data (Thematic Map
For the standard processing data of (per Data), the satellite image on the map (U
As information (system information) for projecting on the TM coordinate system, for example, a geometric transformation formula including a plurality of parameters is included. However, in geometric transformation based on system information (system geometric transformation), since a considerable error occurs, it is generally rarely used. The Landsat TM data that can be obtained from the Space Development Agency is an image obtained by applying a predetermined geometric correction to an actual satellite image and converting the image. However, the geometric correction accuracy of this converted image is not necessarily high. Therefore, conventionally, when geometrical accuracy is required, a precise geometric correction method (including orthographic projection) based on affine transformation using a plurality of GCPs (ground control points) is generally adopted. is there. This GCP
In the precise geometric correction method using (ground control point), the affine transformation formula is statistically estimated by using the least squares method from the coordinates on the map corresponding to the plurality of GCP satellite images visually identified. To do. It is said that it is necessary to acquire as many GCPs as possible in order to statistically increase the precision of the precision geometric correction, because an error is unavoidable in determining the GCP in pixel units. However, there is a problem that the acquisition of GCP is inefficient and high accuracy cannot be expected unless a skilled person performs it.

Therefore, the inventors of the present invention use, as a template (simulation image), a direct irradiance image of the sun, which can be created by using a digital elevation model as a method for automatically determining a large number of GCPs, and after this system geometric transformation. A method for locally matching satellite images (simulation image method) was proposed ("Photogrammetry and remote sensing" Vol. 37, No. 37, entitled "Identification of ground control points of satellite images using digital elevation model"). 6, 1998).

[0004]

As a result of applying the simulation image method previously proposed by the inventor to a large number of Landsat TM images, if this method is adopted, the system geometric transformation is performed when the target area is not so wide. It was found that the size and direction of the positional deviation of the satellite image due to were almost independent of the place. It is speculated that this is because the attitude control of Landsat equipped with a star sensor is performed with extremely high accuracy.

However, in this conventional simulation image method, although the satellite image can be locally projected (matched) on the map (the coordinate system of the satellite image and the coordinate system of the map can be locally matched), the target region The satellite image could not be projected on the map as a whole.

An object of the present invention is to provide a precise geometric correction method for a Landsat TM image and a precise geometric correction for a satellite image capable of performing a precise geometric correction for the entire target area without acquiring a large number of GCPs. To provide a method.

[0007]

In the present invention, the parameters of the geometric transformation formula used in the system geometric transformation are optimized.
At a specific level, the translation parameters are optimized.
At this time, parameters are selected (optimized) using the correlation coefficient between the satellite image and the simulation image as an evaluation function. Since only two parallel movement parameters are specifically changed, it is easy to search for an optimum solution, and fine adjustment within one pixel is possible.

In view of this, the present invention uses a geometric transformation equation including a plurality of parameters to represent the relationship between the coordinates on the map and the coordinates of the Landsat TM image, and the coordinates of the Landsat TM image and the coordinates of the map. An object of the improvement is a precision geometric correction method for a Landsat TM image, which geometrically corrects a Landsat TM image in order to match. In the present invention, a geometric transformation equation including a translation parameter included in the system information provided in association with the Landsat TM image is used as the geometric transformation equation.

For example, the geometric conversion formula contained in this system information is the following formula.

P-p ₀ = [cos θ ₀ * (x-x ₀ )-
sin θ ₀ * (y−y ₀ )] / Δ l−l ₀ = [− sin θ ₀ * (x−x ₀ ) −cos θ ₀
* (Y−y ₀ )] / Δ In the above formula, x ₀ and y ₀ are the center coordinates of the map coordinate system, x and y are the coordinates at any position in the map coordinate system, and p ₀ and l ₀ is the center coordinate (pixel number and line number) of the coordinate system of the Landsat TM image, p and l are the coordinates (pixel number and line number) of any position on the Landsat TM image, and θ ₀ is the map Is the rotation angle of the image coordinate system with respect to the coordinate system of, and Δ is the spatial resolution of the Landsat TM image.

In the present invention, the direct illuminance image of the sun created using the digital elevation model is used as a simulation image, and the appropriate parallel movement is performed by using the correlation coefficient of the pixel values of the simulation image and the Landsat TM image as an evaluation function. Find the parameter (find the parameter that gives the optimum value of the evaluation function). Here, the direct incident illuminance image of the sun is an image including the difference between the direct incident illuminance of the sun obtained from the position information of the sun and the known elevation information. The direct incident illuminance cosβ of the sun is cosβ = cosθ * cos (e) + s
It can be calculated based on the expression in (e) * cos (φ-A). Where θ is the zenith angle, A is the azimuth angle, e is the inclination of the ground surface, and φ is the inclination azimuth. The pixel value is a digital number proportional to the brightness of one pixel in the image, and the pixel value of the Landsat TM image is
It is a value proportional to the brightness of one pixel, and the pixel value of the direct-incidence illuminance image is a value proportional to the illuminance of one pixel.

The correlation coefficient is an index for evaluating the relationship between the respective pixel values in the Landsat TM image and the simulation image, and is a numerical value between -1 and +1. The correlation coefficient can be obtained based on a general formula. Then, in order to use the correlation coefficient as an evaluation function, first, a rough shift amount between the Landsat TM image and the simulation image is obtained by visual inspection or the like, and one image with respect to the other image is centered on the shift amount. The above-mentioned correlation coefficient in a specific pixel area (for example, an area of 100 × 100 pixels at the center of the image) is obtained by moving the pixel direction and the line direction by a predetermined amount. Then, the amount of movement that maximizes this correlation coefficient (it becomes the optimum value) is determined, the amount of deviation between the two images is calculated from this amount of movement, and the value of the parallel movement parameter is set so that this amount of deviation can be corrected. decide.

Directly, p ₀ and l ₀ may be parameters optimized or optimized using an evaluation function.
However, specifically, the above two equations are simplified as follows.

P = ax-by + c l = bx-ay + d where a to d are variables. Then, c and d in this equation are set as parallel movement parameters. And these parameters, the correlation coefficient of the pixel values of the two images as an evaluation function,
The evaluation function is set to have an optimum value, and the above-mentioned p ₀ and l ₀ are obtained.

If the geometric conversion is performed for all the coordinates of the target area according to the above-mentioned geometric conversion equation using the parallel movement parameters determined in this way, it is possible to obtain a high accuracy over the entire target area without acquiring a large number of GCPs. Geometric correction can be performed.

If the value of the correlation coefficient in the pixel direction and the value of the correlation coefficient in the line direction when the parallel movement parameters are changed is the maximum value, the optimum value of the parallel movement parameter is within one pixel. Fine adjustment is possible.

To express the features of the present invention more abstractly,
The feature of the present invention is that the corresponding map information corresponding to the feature on the image among the map-related information obtained in relation to the feature (for example, the above-mentioned brightness) on the image appearing in the entire Landsat TM image and the map information on the map. (For example, the parallel movement parameter is optimized while performing evaluation using an evaluation function that quantitatively evaluates the spatial consistency with the above-mentioned direct incident illuminance of the sun) without being statistically biased. Here, one of the "evaluation functions for quantitatively evaluating spatial consistency without bias" is the pixel value based on the brightness on the TM image and the sun on the direct incident illuminance image (simulation image) of the sun. Is a correlation coefficient with the pixel value based on the direct incident illuminance of. Other such evaluation functions can be used. For example, if the feature on the TM image is the difference between the land and the sea appearing on this image, and the corresponding map information is the difference between the land and the sea shown on the map, then as an evaluation function, Use the identification rate obtained by determining the correctness of the identification data whether the multiple evaluation points randomly selected for each category stratified by the distance from the coastline are land or sea from the corresponding map information. You can It should be noted that "stratification" means that in statistical processing, the population is divided into some subgroups (layers or categories) based on prior information (distance from the coastline). Even when the identification rate is used as the evaluation function, the parameters are determined so that this evaluation function has an optimum value or an appropriate value.

If the features of the present invention are further abstractly (generally) expressed as a precise geometric correction method for satellite images, they can be expressed as follows.

That is, according to the present invention, the coordinate of the satellite image and the coordinate of the map are matched by using the geometric transformation formula constituted by including a plurality of parameters in order to represent the relationship between the coordinate on the map and the coordinate of the satellite image. Targeting the precise geometric correction method of the satellite image that geometrically corrects the satellite image so that the features on the image appearing in the entire satellite image and the map features out of the map-related information obtained in relation to the map information on the map It is characterized in that a plurality of parameters are optimized by using an evaluation function that quantitatively evaluates the spatial consistency with the corresponding map information corresponding to.

Further, in order to facilitate the parameter determination, the parameter may be determined through two steps.
First of all, on the map, the correctness of the judgment as to whether the multiple evaluation points randomly selected for each category stratified by the distance from the coastline on the Landsat TM image (satellite image) is land or sea is displayed on the map. Judgment is made based on the land and sea data shown, and the discrimination rate is obtained. Then, using the identification rate as an evaluation function, appropriate values of a plurality of parameters are obtained. Next, this appropriate value is input to the geometric transformation formula to geometrically correct the satellite image to determine the geometrically corrected satellite image. After that, the direct irradiance image of the sun created using the digital elevation model is used as a simulation image, and the correlation coefficient of each pixel value of the simulation image and the geometrically corrected satellite image is used as an evaluation function. Optimize the value. In this way, it is not necessary to visually check the deviation between the Landsat TM image (satellite image) and the simulation image, and therefore the parameter determination can be easily automated.

[0021]

BEST MODE FOR CARRYING OUT THE INVENTION Embodiments of the present invention will be described in detail below with reference to the drawings. FIG. 1 is a diagram for explaining the flow of data processing of the method according to the first embodiment of this invention. In this embodiment, a Landsat TM image is treated as a satellite image. The information about the map projection ancillary record of the reader file attached to the standard processing data included in the system information of Landsat TM data includes:
The center coordinates on the map that is the reference of the center (center) of the shooting scene, that is, UTM coordinates (x ₀ , y ₀ ), the rotation angle (orientation angle θ ₀ ) of the TM image to be transmitted with respect to the UTM coordinate system, and the spatial resolution (Δ )It is included. Spatial resolution is the smallest unit of object that can be resolved by the sensor that captured the image. If these pieces of information are used, the UTM coordinates x and y on an arbitrary map are converted into positions (lines and pixels) on the satellite image (so-called inverse conversion) using the following geometric conversion formulas (1) and (2). )can do. These equations (1) and (2) are attached to the system information.

P−p ₀ = [cos θ ₀ * (x−x ₀ ) −sin θ ₀ * (y−y ₀ )] / Δ (1) 1−10 = [− sin θ ₀ * (x−x ₀ ) -Cos θ ₀ * (y−y ₀ )] / Δ (2) In the above formula, x ₀ and y ₀ are the center coordinates of the map coordinate system, and x and y are the coordinates of any position of the map coordinate system. Coordinates, p ₀ and l ₀ are the center coordinates (pixel number and line number) of the coordinate system of the Landsat TM image, and p and l are coordinates (pixel number and line number) at any position on the Landsat TM image. Where θ ₀ is the orientation angle and Δ is the spatial resolution.

The above equation can be simplified by the affine transformation equation as follows.

P = ax-by + c l = bx-ay + d where a to d are variables. In this case, c and d are translation parameters. In the system geometric transformation, all the coefficients in the above equation are determined. In the method of the present embodiment, two parameters c and d, which mean parallel movement on image coordinates, are considered as variables. The same is true when p ₀ and l ₀ in the above equations (1) and (2) are considered as parameters (variables).

In this embodiment, the correlation coefficient between a Landsat TM image geometrically corrected through a predetermined geometrical conversion (converted satellite image) and a direct incidence illuminance image created by simulation is used as an evaluation function for geometrical correction. To determine the optimum value of the parameter. Since the correlation coefficient becomes larger as the positioning becomes correct, it becomes possible to optimize the parameters of the geometric transformation formula used for the precise geometric correction.
The optimization is realized by searching the parameters so that the correlation coefficient becomes large.

In this embodiment, since the orthographic projection is included in the geometric transformation using the above formula, the apparent position on the satellite image obtained by adding the positional deviation based on the undulation to the pixel value obtained by the above formula is calculated. Uses data. The geometric transformation used to find the optimal solution shall include orthographic projection. That is, the positional deviation Δp based on the undulations is added to the pixel value p obtained by substituting the map coordinates (x, y) into the equations (1) and (2). However, in order to exactly obtain the positional deviation Δp based on the undulations, it is necessary to perform calculation in consideration of the curvature of the earth, which increases the load on the computer. The calculation time can be shortened by utilizing the fact that the positional deviation due to the strict calculation is about 10% larger than the positional deviation when the curvature of the earth is not considered. The positional deviation due to the undulations of the point directly below the satellite may be estimated from the trajectory of the missing data in the original image. Regarding this point, Vol. Of "Photogrammetry and remote sensing". Three
7, No. 4, page 12-22, entitled "Orthographic Projection of Landsat TM Image and Its Evaluation".

The direct incident illuminance cos β of the sun required to obtain a simulation image is given by the following equation as a function of the sun position (zenith angle θ and azimuth angle A), the inclination e of the ground surface and the inclination azimuth φ. .

Cos β = cos θ cos (e) + sin
θsin (e) cos (φ-A) Here, for the data of the sun position (zenith angle θ and azimuth angle A), the sun position attached to the standard processing data is used.
A numerical elevation model is required to obtain the slope e and the inclination azimuth φ. As the digital elevation model, a numerical map 50 m (elevation) issued by the Geospatial Information Authority of Japan is used after projection conversion (UTM coordinate system) and rearrangement (30 m). The numerical map 50m mesh (elevation) issued by the Geographical Survey Institute is
The grid points (mesh) are composed of evenly spaced latitudes and longitudes, and the size is about 50 m. Therefore, in order to use this for image processing of Landsat TM, projection transformation (UTM coordinate system) and rearrangement (spatial resolution 30 m)
Is needed. For details of this point, see the Remote Sensing Society of Japan, Vol. 21, No. Two
Pp. 150-157 of the article entitled "Evaluation Method of Interpolation Method Used for Projection Transformation of Digital Elevation Model". In the present embodiment, the inclination e and the inclination azimuth φ are obtained from the projection transformation (UTM coordinate system) and rearrangement (spatial resolution 30 m).

In the present embodiment, the target areas are five locations in Iwate prefecture shown in FIG. 2 (areas in which numbers 1 to 5 are written in black circles). A numerical elevation model was created with the size of the satellite image cut out for each point as vertical 430, horizontal 430, and spatial resolution of 30 m.

As the satellite image, Landsat TM data (pass 102, row 32) taken on June 12, 1995 was used. The solar altitude at the time of shooting was 58 degrees, and the azimuth was 112 degrees clockwise from north. Using this information, a simulation image of the direct incident illuminance of the sun was created. Point 1
A simulation image of the above is shown in FIG.

First, a rough alignment was performed by visually comparing the satellite image after system geometric transformation and the simulation image. As a result, a positional shift of 40 pixels in the pixel direction and 27 pixels in the line direction was confirmed. Therefore, the converted Landsat TM image is pixel-wise (one pixel at a time) in pixel and line directions with respect to the simulation image, centering on the positional deviation visually confirmed (40 pixels in the pixel direction, 27 pixels in the line direction). Then, the correlation coefficient between the pixel value of, for example, 100 × 100 pixels located at the center of the Landsat TM image and the pixel value of the corresponding pixel of the simulation image was obtained. The results are shown in Table 1 below.

[0032]

[Table 1] The numbers on the horizontal axis in Table 1 above mean ± 2 pixel shifts with 40 pixels in the pixel direction, and the numbers on the vertical axis mean ± 2 pixel shifts with 27 pixels in the line direction. The numbers in the table are correlation coefficients of pixel values when moving in pixel units. The results in Table 1 show the evaluation function. From Table 1, it is found that the correlation coefficient is 0.6 at a position shifted by 40 pixels in the pixel direction and 27 pixels in the line direction.
It can be seen that 95 is the largest. That is, the optimum value of the evaluation function is 40 pixels in the pixel direction and 2 in the line direction.
It is to shift by 7 pixels. At this stage, the parallel movement parameters c and d described above are obtained so that such a shift amount can be obtained.
Can be determined.

FIG. 4 shows the satellite image (Band 5) after conversion at point 1. It can be confirmed that the shadows of the two are the same as in the simulation image shown in FIG.

Next, a method of determining the optimum value of the parallel movement parameter from the geometric conversion formula having the above-mentioned positional deviation as the origin, which is performed for fine adjustment within one pixel, will be described. Pixel values of 100 × 100 pixels in the central area of a satellite image (Landsat TM image with 40 pixels in the pixel direction and 27 pixels shifted in the line direction) created by changing the parallel movement parameters c and d, and the corresponding simulation The correlation coefficient of the pixel value of the pixel of the image was calculated. Examples of the results are shown in FIGS. 5 and 6. FIG. 5 shows changes in the correlation coefficient when the parallel movement parameter c that determines the movement amount in the pixel direction is changed, and shows the phase when the parallel movement parameter d that determines the movement amount in the line direction is changed. The change in the number of relationships is shown in FIG. Comparing FIGS. 5 and 6, the change is more apparent when the parameter c is changed than when the parameter d is changed, but both have a clean symmetry. By searching two parameters two-dimensionally, the optimum parameter can be determined in the order of one digit smaller than one pixel. It seems that this is because the accuracy of the result is improved because a large number of points can be used for the calculation of the correlation coefficient. That is, the values of the respective parameters c and d (0.2 and 0.2 in this example) when the correlation coefficient becomes the maximum value by changing the parameters c and d can be determined as the respective optimum values. If the image range for calculating the correlation coefficient is set to a range in which the ground cover is uniform and the undulation is relatively large, a correlation coefficient exceeding 0.8 can be obtained. However, even in such a case, a large change is not seen in the information of the positional deviation, and therefore it is preferable to set the range of the correlation coefficient to the entire range in the sense of simply increasing the number of samples.

Table 2 shows the results obtained by thus determining the optimum amount of parallel movement for each of the points 1 to 5 in FIG.

[0036]

[Table 2] In Table 2, c and d on the horizontal axis indicate parameters, and r indicates the value of the correlation coefficient. The selected five points occupy a fairly wide range, but the difference between the optimum translation parameters c and d is 1.
It was about one pixel. In this range, even if c and d are fixed, the number of pixels is less than 1 pixel and less than about 1 pixel in the line direction. This is considered to indicate that the accuracy of geometric correction based on system information is fairly good. In addition, systematic characteristics such as a large deviation toward the minus side in the pixel direction at the west side points (points 2 and 5 in FIG. 1) and a large deviation toward the plus side in the line direction can be seen. This suggests that optimization including rotation is necessary when a larger range is concerned.

In the above description, the optimum values of the parameters c and d have been obtained. However, p ₀ in the above equations (1) and (2)
And l _0, the center of the scene, can be a parameter. In this case, the position shift is adjusted in units smaller than one pixel by changing (δp, δl) the center (p ₀ , l ₀ ) of the scene of the geometric conversion formula. Further, by performing the orthographic projection, it can be expected that an accurate output image is created in a wide area. Furthermore, the calculation of the optimal solution is
This can be done using the Powell method incorporated in IDL (Interactive Data Language) (Ver5.3).

In the example of FIG. 1, the parameters are substituted into the geometric conversion formulas (1) and (2) after the parameters are optimized, and the precise geometric conversion of the Landsat TM image is performed. However, as shown in FIG. 7, when geometrically transforming Landsat TM data, the equations (1) and (2) are used.
Alternatively, the appropriate translation parameter values may be first inserted into these equations to optimize the parameters and, as a result, change the translation parameters.

The present embodiment proposes a simple method (correction system geometric conversion method) for geometrically converting standard processing data of Landsat TM. This modified system geometric transformation method is based on the information about the map projection included in the system information, and considers the transformation that considers the displacement (orthographic projection) due to the parallel movement and the undulation on the ground. Then, the parallel movement parameter is adjusted using the correlation coefficient of the pixel value between the simulation image of the direct incident illuminance of the sun and the satellite image as an evaluation function. Fine adjustment in units of 0.1 pixel is also possible. In the conventional precise geometric correction using GCP, a two-step operation of selecting GCP and statistically estimating the coordinate conversion formula is required, whereas in the method of the present invention, GCP There is an advantage that selection can be omitted.

In the above embodiment, the two parallel movement parameters are optimized, but the parameters that can be optimized are not limited to the parallel movement parameters, and other parameters can be optimized by the same method. Of course, the digitization can be performed.

FIG. 8 is a diagram showing the flow of data processing according to the second embodiment of the present invention. In the above-described first embodiment, the parameters are optimized by visually observing the deviation amount between the satellite image and the simulation image. In this embodiment, instead of visually detecting the amount of deviation,
The amount of global deviation is automatically detected. Therefore, in this embodiment, after detecting the amount of global deviation by a global method (method utilizing the difference between the sea and land), the local method used in the first embodiment described above is detected. The parameters are determined by the method (method utilizing the estimated value of the direct incident illuminance of the sun).

More specifically, in this embodiment, a plurality of evaluation points randomly extracted for each category stratified by the distance from the coastline on the Landsat TM image are land or sea. Whether the judgment is correct or not is judged based on the land and sea data shown on the map, and the discrimination rate is obtained. Then, using this discrimination rate as an evaluation function, an appropriate value of the parallel movement parameter is obtained. Then, an appropriate value is input to the geometric conversion formula to geometrically correct the Landsat TM image to determine a geometrically corrected Landsat TM image. Next, the direct irradiance image of the sun created using the digital elevation model is used as a simulation image, and the correlation coefficient of each pixel value of the simulation image and the geometrically corrected Landsat TM image is used as an evaluation function, and the optimum value of the parallel movement parameter is obtained. To convert.

A global method (method utilizing the difference between the sea and land) will be described below. The most distinguishing feature of satellite images is the difference between land and sea, which can be clearly identified especially in the near infrared band image. The judgment of whether each point on the map coordinates is sea or land can be known by the numerical information issued by the Geographical Survey Institute. To take advantage of this feature, select multiple evaluation points from land and sea areas.
Next, the coordinate of the evaluation point on the satellite image is obtained by using the set geometric correction formula. The applicability is defined as the proportion of evaluation points at which the sea and land that can be determined from the brightness values on these coordinates match the map information. In this method, the setting of evaluation points is important. By including not only the part near the coast but also the part near the coast, it is possible to optimize without depending on the initial value of the parameter. In addition, since the calculation time can be shortened by adjusting the evaluation points, it is effective when targeting a wide area such as the entire scene or when estimating a large number of parameters.

In this embodiment, the satellite image (TM data) from Landsat No. 5 provided by the Space Development Agency of Japan, "photographed on June 12, 1995 (pass 107 / row 3
2) ”was used. The data used here is what is called standard processing data (level 1b). It is (1) adjustment of the sensitivity difference between the detection elements, and (2) projection to the UTM coordinate system using the satellite orbit and attitude information. And rearrangement (system geometric correction) was performed. However, in the system geometric correction, (1) the influence of the undulation of the ground surface is not taken into consideration, (2) the trajectory information is not accurate, (3) the image array does not correspond to the xy axes of the UTM coordinate system ( Due to problems such as orientation angle), it is not possible to accurately overlay the map data.

FIG. 9 shows an image of an example of a visible band image file. The header record (one line) is at the beginning of the file, the first 32 bytes of each record is the record header, and the last 68 bytes is the record trailer.
The information necessary for this analysis is not included. The other portions are image data, but the left and right black portions (luminance value 0) are the portions that have not been photographed that occurred during system correction. In the precise geometric transformation (orthogonal projection), the image is set to the point directly below the satellite at the midpoint of the locus of this boundary for each line.

An information file (read file) is attached to the standard process data in addition to the image file. This file contains the following information necessary for geometric transformation.

Sun altitude θ = 62 degrees Sun direction φ = 119 degrees Orientation angle θ ₀ = 0.1830542 radian Scene center (UTM coordinate system) x ₀ = 534146.1182 m y ₀ = 4463628.9062 m Spatial resolution 28.5 m Also on this image file The center of the scene is p ₀ = 3492.0 pixel l ₀ = 2983.5 line.

The position (p, l) on the image when the UTM coordinates (x, y) are given from these information is expressed by the equations (1) and (2) used in the first embodiment. It is calculated by the following formula.

P = p ₀ + [cos θ ₀ * (x−x ₀ ) −
sin θ ₀ * (y−y ₀ )] / 28.5 1 = l ₀ + [− sin θ ₀ * (x−x ₀ ) −cos θ ₀
* (Y−y ₀ )] / 28.5 Figure 10 shows the satellite data range on the map. In FIG. 10, the area A is the range of the satellite image in FIG. 9, and the area B is
Is the range in which the ground method is implemented, and region C is the range in which the local method is implemented.

The data processing flow shown in FIG. 8 is (A)
It is divided into the creation of evaluation data and (B) estimation of the geometric transformation formula by optimizing the evaluation function. The latter includes a global method using the discrimination between the sea and land and a local method using the direct incident illuminance. Although they can be used independently of each other, in the case of Landsat TM, it is preferable to use a local method by modifying a translation parameter after creating a geometric transformation formula including rotation by a global method. It is rational.

In the data processing flow of FIG. 8, first, evaluation data is created. There are two types of data for evaluation (sea and land discrimination data and direct incident illuminance of the sun), both of which are created from the numerical map 50m mesh (elevation) issued by the Geographical Survey Institute.

FIG. 11 is a cutout of the digital elevation model of a rectangular area B where the range of the target satellite image (for path 107 / row 32) and most of the coastline shown in FIG. 11 overlap.
It is shown in. Using this image, a distance image from the coastline (level 1 where the sea and land are in contact) was created. The area of about 750 m from the coastline is stratified into 15 levels by distance, and 200 pixels (half the sea and half the land) are randomly extracted for each level (category), and a total of 3000 points are evaluated. The data as shown in Table 3 was created as the points.
FIG. 12 shows a part of this. In FIG. 12, the + mark indicates the evaluation point. Further, the data used in the determination of the geometric conversion formula converts the position data indicated by the latitude and longitude into the UTM coordinate system.

[0053]

[Table 3] FIG. 13 shows a direct incident illuminance image (simulation image) of the sun in the area C shown in FIG. 10. This image shows a digital elevation model (UTM coordinate system / spatial resolution 30
m) is created by projection transformation / rearrangement (colinear interpolation) of a 50 m mesh (elevation) of a digital map. The size of the image is 430 × 430. Direct incident illuminance of the sun co
The method of obtaining sβ is the same as in the case of the first embodiment.

Next, the estimation of the geometric conversion formula by optimizing the evaluation function will be described.

A global method is applied to the satellite image of band 4 shown in FIG. The position of the evaluation point on the satellite image is calculated based on the system information, and the case where the brightness value (pixel value) of each evaluation point is less than or equal to the threshold value γ (set 20) is set as the sea area, and the brightness value is greater than the threshold value. Judged to be land area. Then, whether the judgment result (identification data) is correct or not is judged by referring to the map-related information obtained from the above-mentioned numerical map 50 m mesh (elevation), and the discrimination rate is obtained. As a result, 718 points were erroneously determined from 2878 points excluding 122 points outside the range of the satellite image from the evaluation points.

In this example, a wrong number is used as the discrimination rate, and this is used as an evaluation function to increase the evaluation value (to reduce the discrimination rate).
Optimization with respect to Δl) and threshold (γ) reduced the number of errors to 51. Here, the optimum value is Δ
p, Δl, and the threshold value γ were changed according to a known optimization method, and the brightness value (pixel value) and the threshold value were compared with each other. The parameter values at this time are as follows.

Case A: Δp = 40.90, Δl = 27.18,
γ = 28.39 Case A has only parallel movement, but when the optimization was performed to include rotation, the deviation Δθ of the orientation angle was further included, and 40 errors were found. Figure 1
The wrong place and the out-of-range point are shown on the map in 4. What is visible in the north is the evaluation points outside the range of the satellite image. The parameter values at this time are as follows.

Case B: Δp = 40.25, Δl = 27.18,
When the optimization was performed by changing the spatial resolution of 28.5 m, Δθ = 0.000217, γ = 29.08, no error phenomenon was observed. Although we did not check the wrong points individually, there is no systematic tendency, except that there are many lakes and marshes around Mutsu Ogawara. It may be due to the recovery of small islands, estuaries, or ports.

Images were rearranged using the geometric transformation formula obtained by the above-mentioned global method. That is, a rough shift amount between the geometrically corrected Landsat TM image and the simulation image was determined. The rearrangement of this image corresponds to the visual determination of the amount of shift of the center of the scene in the first embodiment. After that, similarly to the first embodiment, the Landsat TM image geometrically corrected in pixel units is moved in the pixel direction and the line direction so that the pixel values of the converted satellite image and the simulation image have different pixel values. The change in the number of relationships was sought. Tables 4 and 5 show the changes in the correlation coefficient obtained when the parameters of the above-mentioned cases A and B were used, respectively. It should be noted that in this example, in order to avoid an error due to ups and downs of the terrain, orthoscopic projection is taken into consideration for the conversion. A co-linear interpolation method was used for the interpolation. For the analysis, the correlation coefficient obtained by shifting the converted satellite image in the pixel direction and the line direction was used. If it is a correct conversion, it takes the maximum value where there is no deviation.

[0060]

[Table 4]

[Table 5] In Tables 4 and 5, the horizontal axis indicates the amount of movement (pixel number) in the pixel direction, and the vertical axis indicates the amount of movement (pixel number) in the line direction. The numbers on the horizontal and vertical axes are the center shifts Δp of the scenes obtained in Case A and Case B.
It means that the image is shifted in the pixel direction and the line direction on a pixel-by-pixel basis based on the values of Δl and Δl (set to 0).

In both case A (only parallel movement) and case B (including rotation), there is no deviation in the line direction, and a maximum deviation of about 2 pixels is observed in the pixel direction. Comparing the two, the positional deviation of the latter is slightly smaller.

Next, the parameters were optimized for each case using the correlation coefficient as an evaluation function. The results are shown below.

Case A: Δp = 38.79, Δl = 27.45,
Correlation coefficient 0.692 Case B: Δp = 39.19, Δl = 27.21, Δθ = −0.00
[0213] Correlation coefficient 0.692 When optimization is performed, the effect of including rotation hardly appears. This is because in the range of the target image (about 10 km), the influence of rotation appears only for several meters (10000 × Δθ). Therefore, a global method is indispensable for obtaining the rotational deviation. Also, when performing optimization by a local method, it may be difficult to find a solution unless the initial value is set appropriately. This problem can be solved by using the parameters obtained by the global method as the initial values.

[0064]

According to the present invention, since an appropriate translation parameter can be obtained by using the correlation function as an evaluation function, geometric transformation can be performed for all coordinates of the target area using this translation parameter. There is an advantage that geometric correction can be performed with high accuracy over the entire target region without acquiring a large number of GCPs.

[Brief description of drawings]

FIG. 1 is a diagram for explaining a flow of data processing of a method according to a first embodiment of this invention.

FIG. 2 is a diagram showing a target area.

FIG. 3 is a diagram showing a simulation image of a point 1 shown in FIG.

4 is a diagram showing a satellite image (band 5) after conversion of point 1 shown in FIG.

FIG. 5 is a diagram showing changes in a correlation coefficient with respect to a parameter c.

FIG. 6 is a diagram showing changes in a correlation coefficient with respect to a parameter d.

FIG. 7 is a diagram showing another data processing flow according to the first embodiment.

FIG. 8 is a diagram showing a flow of data processing according to the second embodiment of this invention.

FIG. 9 is a diagram showing an image of a visible band image file.

FIG. 10 is a diagram showing a target area.

11 is a diagram showing a digital elevation model of a region B in FIG.

FIG. 12 is a diagram showing an example of evaluation points randomly selected.

13 is a diagram showing a direct-incidence illuminance image (simulation image) of area C shown in FIG.

FIG. 14 is a diagram showing evaluation points in which land and sea are erroneously determined.

[Explanation of symbols]

Area A to C

   ─────────────────────────────────────────────────── ─── Continued front page    F-term (reference) 5B050 AA01 BA02 BA07 BA08 BA17                       DA07 EA13                 5B057 AA13 AA14 CA08 CA13 CA16                       CB08 CB13 CB16 CD11 CH01                       DB03 DB09 DC30                 5J070 AC01 AE07 AF08 AH04 AH19                       AJ02 AK22 AK37 BG27                 5L096 AA06 AA09 BA20 DA01 EA15                       FA34 FA66 FA69 FA78

Claims (9)

[Claims] 1. The coordinates of the Landsat TM image and the coordinates of the map are calculated by using a geometric transformation equation including a plurality of parameters to represent the relationship between the coordinates on the map and the coordinates of the Landsat TM image. What is claimed is: 1. A method of precisely correcting a Landsat TM image for geometrical correction to match the Landsat TM image, the method including a translation parameter included in system information provided accompanying the Landsat TM image as the geometric conversion formula. Using the geometric transformation formula, the direct irradiance image of the sun created by using the digital elevation model is used as a simulation image, and the correlation coefficient of each pixel value of the simulation image and the Landsat TM image is used as an evaluation function. Precise geometrical compensation of Landsat TM image, characterized by obtaining the parallel movement parameter Method. 2. A value that maximizes the values of the correlation coefficient in the pixel direction and the correlation coefficient in the line direction when each of the translation parameters is changed is the optimum value of the translation parameter. Claim 1 characterized by the above.
The method for precise geometrical correction of a Landsat TM image according to 1. 3. The coordinates of the Landsat TM image and the coordinates of the map are calculated by using a geometric transformation equation including a plurality of parameters to represent the relationship between the coordinates on the map and the coordinates of the Landsat TM image. What is claimed is: 1. A method of precisely correcting a Landsat TM image for geometrical correction to match the Landsat TM image, the method including a translation parameter included in system information provided accompanying the Landsat TM image as the geometric conversion formula. By using a geometric transformation formula, the spatial relationship between the feature on the image appearing in the entire Landsat TM image and the corresponding map information corresponding to the feature on the image among the map-related information obtained in association with the map information on the map Using an evaluation function for quantitatively evaluating the consistency without statistically biasing, the parallel movement parameters are optimized. Precision geometric correction method of that Landsat TM image. 4. The feature on the image is a pixel value based on luminance, and the corresponding map information is a pixel value based on direct incident illuminance of the sun obtained from position information of the sun and known elevation information, 4. The Landsat T according to claim 3, wherein a correlation coefficient between the pixel value based on the brightness and the pixel value based on the direct incident illuminance of the sun is used as a function.
Method for precise geometric correction of M image. 5. The feature on the image is the difference between the land and the sea appearing on the image, the corresponding map information is the difference between the land and the sea shown on the map, and the evaluation function is , The correctness of the identification data whether a plurality of evaluation points randomly extracted for each category stratified by the distance from the coastline are land or sea, the identification rate obtained by determining by the corresponding map information The method for precise geometry correction of a Landsat TM image according to claim 3, which is used. 6. The geometric transformation formula contained in the system information is p−p ₀ = [cos θ ₀ * (x−x ₀ ) −sin θ ₀ *
(Y−y ₀ )] / Δ 1−l ₀ = [− sin θ ₀ * (x−x ₀ ) −cos θ ₀
* (Y−y ₀ )] / Δ, where x ₀ and y ₀ are the center coordinates of the map coordinate system, and x and y are the coordinates of any position in the map coordinate system. Yes, p ₀ and l ₀ are the center coordinates of the coordinate system of the Landsat TM image, p and l are the coordinates of any position on the Landsat TM image,
θ ₀ is the rotation angle of the coordinate system of the image with respect to the coordinate system of the map, Δ is the spatial resolution of the Landsat TM image, and p = ax−by + c and l = bx are given by the above equations.
The method for precise geometric correction of a Landsat TM image according to claim 1 or 3, wherein the translation parameters are the c and d when simplified to -ay + d (where a to d are variables). 7. The geometric conversion formula contained in the system information is p−p ₀ = [cos θ ₀ * (x−x ₀ ) −sin θ ₀ *
(Y−y ₀ )] / Δ 1−l ₀ = [− sin θ ₀ * (x−x ₀ ) −cos θ ₀
* (Y−y ₀ )] / Δ, where x ₀ and y ₀ are the center coordinates of the map coordinate system, and x and y are the coordinates of any position in the map coordinate system. Yes, p ₀ and l ₀ are the center coordinates of the coordinate system of the Landsat TM image, p and l are the coordinates of any position on the Landsat TM image,
θ ₀ is the rotation angle of the image coordinate system with respect to the map coordinate system, Δ is the spatial resolution of the Landsat TM image, and p ₀ and l ₀ are parameters optimized using the evaluation function. The method for precise geometric correction of a Landsat TM image according to claim 1 or 3, wherein: 8. The coordinates of the satellite image and the coordinates of the map are matched by using a geometric transformation formula configured to include a plurality of parameters to represent the relationship between the coordinates on the map and the coordinates of the satellite image. A satellite image precision geometric correction method for geometrically correcting the satellite image, wherein the image out of the map-related information obtained in relation to the features on the image appearing in the entire satellite image and the map information on the map Precise geometry of satellite images characterized by performing optimization of the plurality of parameters using an evaluation function that quantitatively evaluates spatial consistency with corresponding map information corresponding to the above features, without statistically biasing Correction method. 9. The coordinates of the satellite image and the coordinates of the map are matched by using a geometric transformation formula configured to include a plurality of parameters to represent the relationship between the coordinates on the map and the coordinates of the satellite image. For this purpose, there is provided a satellite image precision geometric correction method for geometrically correcting the satellite image, wherein a plurality of evaluation points randomly extracted for each category stratified by the distance from the coastline on the satellite image is land. Whether or not it is the sea is determined based on the land and sea data shown on the map, the identification rate is obtained, and the appropriate value of the plurality of parameters is obtained using the identification rate as an evaluation function. , Inputting the appropriate value into the geometric conversion formula to geometrically correct the satellite image to determine a geometrically corrected satellite image, and use a direct elevation irradiance image of the sun created using a numerical elevation model as a simulation image, The above The correlation coefficient for each pixel value of the simulation image and the geometric correction satellite images as an evaluation function,
A method for precisely geometrically correcting a satellite image, which comprises optimizing the plurality of parameters.