Next Article in Journal
Review of Automatic Feature Extraction from High-Resolution Optical Sensor Data for UAV-Based Cadastral Mapping
Next Article in Special Issue
Methodology for Detection and Interpretation of Ground Motion Areas with the A-DInSAR Time Series Analysis
Previous Article in Journal
Built-up Area Extraction from PolSAR Imagery with Model-Based Decomposition and Polarimetric Coherence
Previous Article in Special Issue
Ground Subsidence in the Beijing-Tianjin-Hebei Region from 1992 to 2014 Revealed by Multiple SAR Stacks
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Second-Order Polynomial Equation-Based Block Adjustment for Orthorectification of DISP Imagery

1
Guangxi Key Laboratory for Geospatial Informatics, Guilin University of Technology, No. 12, Jian’gan Road, Qixing District, Guilin 541004, China
2
Guangxi Institute of Surveying and Remote Sensing Information, No. 5, Jianzheng Road, Qingxiu District, Nanning 530023, China
3
The Center for Remote Sensing, Tianjin University, No. 92, Weijin Road, Nankai District, Tianjin 300072, China
*
Authors to whom correspondence should be addressed.
Submission received: 4 July 2016 / Revised: 12 August 2016 / Accepted: 17 August 2016 / Published: 22 August 2016
(This article belongs to the Special Issue Earth Observations for Geohazards)

Abstract

:
Due to the lack of ground control points (GCPs) and parameters of satellite orbits, as well as the interior and exterior orientation parameters of cameras in historical declassified intelligence satellite photography (DISP) imagery, a second order polynomial equation-based block adjustment model is proposed for orthorectification of DISP imagery. With the proposed model, 355 DISP images from four missions and five orbits are orthorectified, with an approximate accuracy of 2.0–3.0 m. The 355 orthorectified images are assembled into a seamless, full-coverage mosaic image map of the karst area of Guangxi, China. The accuracy of the mosaicked image map is within 2.0–4.0 m when compared to 78 checkpoints measured by Real–Time Kinematic (RTK) GPS surveys. The assembled image map will be delivered to the Guangxi Geological Library and released to the public domain and the research community.

Graphical Abstract

1. Introduction

Rocky karstification in karst areas (also called karst rocky desertification (KRD)) is considered one of the major factors that contribute to the global carbon balance as a global CO2 sink [1,2,3]. With the increasing interest in global carbon emissions, studies and analyses have compared historical data with current data to discover how rocky karstification contributes to long-term environmental changes over decadal spans.
Guangxi is located in the southwestern karst area in China, and the KRD area is approximately 23,790.80 km2, accounting for 19.8% of the total KRD area in China in 2005. It shrunk to 19,260.00 km2 in 2011, accounting for 16.0% of the total KRD area in China. Although many researchers have investigated the Guangxi KRD area associated with its environmental evolution in recent decades, there have been no investigations or analyses of the KRD area that focused on the early 1960s. Fortunately, declassified intelligence satellite photography (DISP) released to the public domain in February 1995 has provided researchers with a unique opportunity to investigate the KRD in Guangxi in the 1960s. The DISP was collected by the first generation of United States photoreconnaissance satellites between 1960 and 1972 through the systems named CORONA, ARGON, and LANYARD. More than 860,000 images of the Earth’s surface were declassified with the issuance of this executive order and were contracted to the USGS for sale.
However, further processing and application of DISP has resulted in various problems:
(1) The USGS does not provide Chinese users with the parameters required to further process DISP. These parameters include satellite orbit parameters (e.g., inclination, flight height, descent time, etc.) and the camera’s interior orientation parameters (IOP) (e.g., focal length, principal point coordinates, fiducial marks, etc.). This implies that traditional bundle block adjustment based on the photogrammetric collinearity equation is not applicable [4,5].
(2) It is very difficult to obtain sufficient ground control points (GCPs) in the historical DISP imagery due to the time intervals of several decades and cloudy coverage in Southern China. Thus, it is almost impossible to rectify each DISP image on a frame-by-frame base.
For the two reasons above, this paper presents a second order polynomial equation-based rectification model for orthorectification of DISP images. The previous relevant studies on this topic are as follows: Kim et al. utilized a collinearity equation to rectify ARGON imagery from 1963 to study the seasonal variations of glaciers on the Queen Maud Land coast of Antarctica [6]. Zhou and Jezek proposed a collinearity equation-based self-calibration block bundle adjustment method that integrates the bundle adjustment method and satellite orbital parameters, solving interior orientation (including lens distortion) and exterior orientation parameters (EOPs) simultaneously to rectify ARGON images from 1962 and 1963 [4]. The rectified ARGON imagery was employed to mosaic Greenland ice sheets from the 1960s, which were then quantitatively compared to the ice sheet extent over a 30-year interval [5]. Kim and Jezek applied a state-of-the-art digital imaging technology based on an extended block adjustment to rectify ARGON imagery from 1963 that covered Antarctica [7]. They assembled all images into a quality mosaic of coastal Antarctica to study glaciers. In addition, due to the imaging model limitations of high-resolution satellites, such as IKONOS, rational polynomial-based block adjustment, also called rational polynomial coefficient (RPC), was proposed by multiple authors. For example, Tao et al. analyzed the accuracy of orthorectification of a Systeme Probatoire d’Observation de la Terre (SPOT) image and an aerial image using the RPC model [8,9]. Yang suggested that the RPC model can replace the rigorous sensor model orthorectification of SPOT images [10]. Liu developed a stereotaxic method of IKONOS images based on the RPC model [11]. Huang proposed a rational polynomial-based block adjustment for orthorectification of Synthetic Aperture Radar (SAR) images [12]. Grodecki and Gene Dial rectified IKONOS satellite imagery using the RPC method. The RPC model incorporates a priori constraints into the images described by RPC, and multiple independent images can be added in accordance with the needs of users [13]. However, the RPC model requires a number of GCPs, and the computation is very time consuming. Therefore, the RPC method is not applicable to DISP images that the USGS provides because the imaging model of DISP was not provided by the USGS. Additionally, few GCPs are available in the study area. Thus, this paper presents an effective and simple mathematical model for geometric rectification of DISP images, considerably improving the computational effectiveness.

2. The Second Order Polynomial Equation-Based Rectification Model Method

2.1. Polynomial Equation-Based Block Adjustment Model

The objective of polynomial equation-based block adjustment is to tie overlapping images together without the absolute need for ground control points in each image and obtain coordinates of tie points and conversion parameters for rectification. Since the study area is a karst landform with a large wavy terrain and large elevation differences, the relief displacement is large. For correction of relief displacement, relief displacement is introduced into the block adjustment model shown below.
Figure 1 shows the imaging geometry of DISP from the CORONA mission. S W V U ( W V U ) is a camera coordinate system, o x y is an image plane system, and O X Y Z is a geographic coordinate system. There is relief displacement ( Δ h ) in the imaging process; therefore, the relief displacement must be corrected in the rectification process. First, distortion caused by elevation differences should be corrected. Then, other distortions should be corrected by utilizing a polynomial model, and the reverse is true in the resampling process.
Since the relief displacement only occurs in the direction of scanning, CORONA images are panoramic camera images, and the panoramic projection scan direction is the x -direction. Therefore, as shown in the imaging equation above, there is no relief displacement in the y -direction. The relief displacement correction functions are as follows:
Δ h = Z h / M
{ Δ x = x Z / M Δ y = y Z / M = 0
where x and y are image coordinates; Δ x and Δ y are image distortions in the x - and y -directions, respectively, caused by elevation differences; Z is the elevation; h is the distance from the image point to the nadir point; and M is the satellite flight altitude. Since the relief displacement occurs in the direction of scanning and the KH-4A/B’s images are panoramic camera images, the images can be rectified using the second-order polynomial equation-based model.

2.1.1. Traditional Second-Order Polynomial Equation

The traditional second order polynomial model has been widely applied for image rectification. This paper extends the traditional equation into a block situation by adding tie points, which tie overlapping images together. With the extended model, the 2D coordinates of tie points and the coefficients of second order polynomial equations are solved. Furthermore, these parameters are used for orthorectification of DISP imagery without the absolute requirement of at least six GCPs in each DISP image.
The traditional second order polynomial equations are expressed as follows [14]:
{ x + Δ x = a 0 + a 1 X + a 2 Y + a 3 X Y + a 4 X 2 + a 5 Y 2 y + Δ y = b 0 + b 1 X + b 2 Y + b 3 X Y + b 4 X 2 + b 5 Y 2
where α = ( a 0 , a 1 , a 2 , a 3 , a 4 , a 5 ) T and β = ( b 0 , b 1 , b 2 , b 3 , b 4 , b 5 ) T are coefficients; x and y are image coordinates; Δ x and Δ y are image distortions in the x - and y -directions, respectively; and X and Y are 2D coordinates in a given map coordinate system.
For a given GCP, Equation (3) can be linearized using a Taylor series and is expressed as follows:
{ v x = Δ a 0 + X Δ a 1 + Y Δ a 2 + X Y Δ a 3 + X 2 Δ a 4 + Y 2 Δ a 5 l x v y = Δ b 0 + X Δ b 1 + Y Δ b 2 + X Y Δ b 3 + X 2 Δ b 4 + Y 2 Δ b 5 l y
where Δ a i ( i = 0 , 1 , , 5 ) and Δ b i ( i = 0 , 1 , , 5 ) are correction terms of coefficients; v x ,   v y are residuals; X and Y are 2D coordinates of GCPs; and l x and l y are constants expressed by Equation (5).
{ l x = x ( a 0 + a 1 X + a 2 Y + a 3 X Y + a 4 X 2 + a 5 Y 2 ) l y = y ( b 0 + b 1 X + b 2 Y + b 3 X Y + b 4 X 2 + b 5 Y 2 )
As shown in Equation (4), one GCP only establishes two observations, but Equation (4) has 12 unknown parameters. Therefore, six GCPs, which establish 12 observation equations, are needed to solve the 12 coefficients that are used for to rectify a single image. Generally, more than six GCPs are observed in each image to establish more than 12 observation equations. The least-squares estimation is employed to calculate the 12 coefficients. Mathematically, the solution can be described as follows.
Assuming that N GCPs (N ≥ 6) are observed, the observation equations are expressed in matrix form as follows:
V = A α L
where:
V = [ V x 1 V y 1 ... ... V x N V y N ] T ,
A = ( 1 X 1 Y 1 X 1 Y 1 X 1 2 Y 1 2 0 0 0 0 0 0 0 0 0 0 0 0 1 X 1 Y 1 X 1 Y 1 X 1 2 Y 1 2 1 X N Y N X N Y N X N 2 Y N 2 0 0 0 0 0 0 0 0 0 0 0 0 1 X N Y N X N Y N X N 2 Y N 2 ) ,
α = ( a 0 a 1 a 2 a 3 a 4 a 5 b 0 b 1 b 2 b 3 b 4 b 5 ) T
and:
L = [ l x 1 l y 1 ... ... l x N l y N ] T .
The least-squares estimation, i.e., V T P V = min , gives the solutions of the coefficients of the second order polynomial equation below:
α = ( A T A ) 1 A T L
We can further obtain the following expressions from Equation (7):
a i = a i 0 + j = 1 N i t e Δ a i j ( i = 1 , , 5 ;   j = 1 , , N i t e )
b i = b i 0 + j = 1 N i t e Δ b i j ( i = 1 , , 5 ;   j = 1 , , N i t e )
where a i 0 , b i 0 are initial values; Δ a i j , Δ b i j are increases during each iteration; and N i t e is the number of iterations.

2.1.2. The Second-Order Polynomial Equation-Based Rectification Model

As mentioned above, due to the shortage of GCPs in each of the DISP images, the tie points (TPs) must be identified to tie images with the same overlapping areas. Under this condition, the TPs whose XY-coordinates are unknown are introduced into the traditional second order polynomial equation. This extended model is called the second-order polynomial equation-based rectification model (2OPE-RM) in this paper (see Figure 1). Equation (3) is extended with considering TPs as unknown parameters and linearized into the following form:
{ v x = Δ a 0 + X Δ a 1 + Y Δ a 2 + X Y Δ a 3 + X 2 Δ a 4 + Y 2 Δ a 5 + ( a 1 + a 3 Y + 2 a 4 X ) Δ X + ( a 2 + a 3 X + 2 a 5 Y ) Δ Y l x v y = Δ b 0 + X Δ b 1 + Y Δ b 2 + X Y Δ b 3 + X 2 Δ b 4 + Y 2 Δ b 5 + ( b 1 + b 3 Y + 2 b 4 X ) Δ X + ( b 2 + b 3 X + 2 b 5 Y ) Δ Y l y
Then, Equation (10) can be rewritten as follows:
{ v x = Δ a 0 + X Δ a 1 + Y Δ a 2 + X Y Δ a 3 + X 2 Δ a 4 + Y 2 Δ a 5 + f 1 Δ X + f 2 Δ Y l x v y = Δ b 0 + X Δ b 1 + Y Δ b 2 + X Y Δ b 3 + X 2 Δ b 4 + Y 2 Δ b 5 + g 1 Δ X + g 2 Δ Y l y
where:
f 1 = a 1 + a 3 Y + 2 a 4 X ,
f 2 = a 2 + a 3 X + 2 a 5 Y ,
g 1 = b 1 + b 3 Y + 2 b 4 X ,  and
g 2 = b 2 + b 3 X + 2 b 5 Y .
The symbols above are the same as those in Equation (10).
Additionally, assuming that there are N GCPs (N ≥ 6), M TPs in t images are collected at the GCPs. Similarly, Equation (10) can be expressed in matrix form as follows:
V = A α + B β L
where:
V = [ v x 1 G C P v y 1 G C P ... ... v x n G C P v x n G C P | v x 1 T P v y 1 T P ... ... v x M T P v x M T P ] T ,
α = ( Δ α 0 1 ,   Δ α 1 1 ,   Δ α 2 1 ,   Δ α 3 1 ,   Δ α 4 1 ,   Δ α 5 1 ,   ,   Δ α 0 t ,   Δ α 1 t ,   Δ α 2 t ,   Δ α 3 t ,   Δ α 4 t ,   Δ α 5 t ) ,
β = ( Δ X 1 1 ,   Δ Y 1 1 , , , Δ X M t ,   Δ Y M t ) ,
A = ( 1 X 1 1 Y 1 1 X 1 1 Y 1 1 X 1 1 2 Y 1 1 2 0 0 0 0 0 0 0 0 0 0 0 0 1 X N 1 i Y N 1 i X N 1 i Y N 1 i X N 1 1 2 Y N 1 1 2 } The  1 st  image with  N 1  GCPs 1 X 1 i Y 1 i X 1 i Y 1 i X 1 i 2 Y 1 i 2 0 0 0 0 0 0 0 0 0 0 0 0 1 X N i i Y N i i X N i i Y N i i X N i i 2 Y N i i 2 } The  i th  image with  N i  GCPs 1 X 1 N Y 1 N X 1 N Y 1 N X 1 N 2 Y 1 N 2 0 0 0 0 0 0 0 0 0 0 0 0 1 X N n N Y N n N X N n N Y N n N X N n N 2 Y N n N 2 } The  N th  image with  N n  GCPs )
and:
B = ( f 1 1 f 2 1 f 2 M 1 f 2 M 1 } The  1 st  image with  M 1  TPs f 1 i f 2 i f 1 M i f 2 M i } The  i th  image with  M i  TPs f 1 M f 2 M f 1 M n f 2 M n } The  M th  image with  M n  TPs )
Equation (12) is the 2OPE-RM model derived in this paper. Relative to the traditional model in Equation (6), this model introduces TPs as unknown parameters.
Equation (12) is usually solved using least-squares estimation, which is expressed as follows:
Φ = V T V = min
With least-squares estimation, the normal equation matrix can be written as follows:
( A T A A T B B T A B T B ) ( δ α δ β ) = ( A T L B T A )
Thus, the solution of the unknown parameters is given by Equation (15):
{ Δ a = ( Q α α A T L + Q α β B T L ) Δ b = ( Q β α A T L + Q β β B T L )
where Q i j ( i , j = 1 , 2 ) gives the components of the covariance matrix, which is the inverse of the normal matrix, as shown in Equation (16):
Q i j = ( A T A A T B B T A B T B ) 1 = ( Q α α Q α β Q β α Q β β ) ( i , j = 1 , 2 )
The coefficients of the 2OPE-RM in each image and the 2D coordinates (XY) of each TP are as follows:
a i = a i 0 + j = 1 N i t e Δ a i j ( i = 1 , , 5 ;   j = 1 , , N i t e )
b i = b i 0 + j = 1 N i t e Δ b i j ( i = 1 , , 5 ;   j = 1 , , N i t e )
X i t i = X i 0 + j = 1 N i t e Δ X i j ( i = 1 , , 5 ;   j = 1 , , N i t e ; t i = 1 , , t )
Y i t i = Y i 0 + j = 1 N i t e Δ Y i j ( i = 1 , , 5 ;   j = 1 , , N i t e ; t i = 1 , , t )
where X i ,   Y i are coordinates of the i-th TP in image ti; Δ X i ,   Δ Y i are increases in X i and Y i ; a i 0 and b i 0 are initial values; and Δ a i and Δ b i are increases in the coefficients in each iteration.
As shown in Equation (12), each image has 12 unknown parameters ( a i , b i ; i = 0 , 1 , , 5 ) , and each TP has two unknown parameters (XY-coordinates). Two equations can be established for each GCP or TP. Moreover, the TPs and/or GCPs should be well distributed in each image. For example, there are four images, 12 GCPs, and nine TPs in Figure 2. The four images imply that there are 48 unknown parameters. The 12 GCPs can be used to establish 42 observation equations (i.e., seven GCPs in Image #1 can be used to establish 14 observations, three GCPs in Image #2 can be used to establish six observations, six GCPs in Image #3 can be used to establish 12 observations, and five GCPs in Image #4 can be used to establish 10 observations). The nine TPs can be used to establish 34 observation equations (i.e., three TPs in Image #1 can be used to establish six observations, three TPs in Image #2 can be used to establish six observations, six TPs in Image #3 can be used to establish 12 observations, and 5 GCPs in Image #4 can be used to establish 10 observations). With this model, we have 76 (76 = 42 + 34) observations and 66 (66 = 48 + 18) unknown parameters. Thus, 2OPE-RM does not require each DISP image to have more than six GCPs.
The accuracy of the adjustment computation is evaluated using Equation (21):
δ o = V T V r
where δ o is the standard deviation of the unit weight, V is the matrix of residuals, and r is the number of redundant observations. Thus, the standard deviations of individual unknown parameters can be calculated as follows:
δ X i = δ o Q X i
To evaluate the accuracies of TPs, assuming that there are n TPs, the average of δ X i is as follows:
μ X = 1 n δ X i
where n is the number of TPs.

2.2. Orthorectification of DISP Images

With the established model and the coefficients determined in Section 2.1, each original DISP image can be orthorectified. The steps are as follows:
1.
Step 1: Determination of the Rectified Image Size
To properly establish the storage space of the orthorectified image, the size of the resulting image (upper left, lower left, upper right, and lower right) must be determined in advance. This procedure is proposed as follows.
  • Determination of the four corner coordinates: The four corner coordinates of the original image are projected into the UTM coordinate system. Then, eight coordinates are obtained:
    ( X u l , Y u l ) , ( X l l , Y l l ) , ( X u r , Y u r ) , ( X l r , Y l r ) .
The maximum and minimum values of X and Y ( X min , X max , Y min , and Y max ) are calculated from the eight coordinates above to constitute four coordinate pairs. These pairs are the map coordinates of the four boundaries of the resulting image’s scope.
X min = min ( X u l , X l l , X u r , X l r ) , X max = min ( X u l , X l l , X u r , X l r ) Y min = min ( Y u l , Y l l , Y u r , Y l r ) , Y max = min ( Y u l , Y l l , Y u r , Y l r )
  • Determination of the resulting image’s size: The size of the resulting image can be determined by M and N as follows:
    M = Y max Y G S D + 1 , N = X min X G S D + 1
    where M = r o w , N = c o l , and Y G S D , X G S D are the ground-sampled distances (GSD) in the resulting image.
2.
Step 2: Coordinate Transformation
Because the orthorectification model only expresses the relationship between the original coordinates ( x o r i , y o r i ) and ground coordinates ( X g r o , Y g r o ) , the ground coordinates should be transformed into the coordinates of the resulting image ( x r e , y r e ) as follows:
x r e = Y max Y g r o Y G S D + 1 ; y r e = X g r o X min X G S D + 1
where Y g r o , X g r o are the ground coordinates of the pixel after rectification.
3.
Step 3: Orthorectification
The calculation of the geographic coordinates of individual pixels, resampling of the original image, and registration of the chosen map coordinates system are carried out as follows:
  • The process can be applied to any point P ( I , J ) in the resulting image with image coordinates ( I , J ) .
  • In accordance with image coordinates ( I , J ) and GSD, calculate the geographic coordinates ( X , Y ) .
  • Compute the image coordinates ( i , j ) of point P in the original image using Equation (5).
  • Calculate the gray value g o r i via bilinear resampling interpolation.
  • Assign the gray value g o r i to point P as g r e s in the resulting (rectified) image/pixel.
The above procedure is then repeated for each pixel that must be rectified until the entire image is completely rectified.

2.3. Data Set

  • Study area
The study area is located in Guangxi, China, spanning from 20.54°N to 26.24°N latitudes and 104.26°E to 112.04°E longitudes (Figure 3) and encompassing 23,790.8 km2. The study area is in the south central subtropics of China.
  • Dataset
DISP imagery: In total, 444 DISP images from five orbits of different missions, including the CORONA 1035-1 Mission (24 images) on 25 September 1966, the CORONA 1102-2 Mission (48 images) on 18 December 1967, and the 1106-1/2 Mission (39 images) on 7 February 1969, were purchased from the USGS (Figure 4).
Aerial photos: Five aerial photos with film formats of 18 × 18 cm2 from 1961 were acquired at a photographic scale of 1:14,000. Each photo covers approximately 6.35 km2. Five aerial photos were purchased from the Guangxi Bureau of Geospatial Information, China.
Coordinate data of GCPs: The coordinate data associated with GCPs in the KRD area were collected from Google Earth.

3. Results and Accuracy Analysis

3.1. Image Preprocessing

The DISP film was scanned into digital images, producing film-grain noise and resulting in image quality degradation. Many noise filters have been used in the public domain. However, most of these approaches are either time consuming, because of complex modelling, or they erroneously remove geophysical features because of noise in the overall image. The filter algorithm developed by Zhou et al. was used to remove noise in this study [5]. One of the advantages of the algorithm is that it avoids the problems noted above because this approach performs statistical calculations within variable-size and variable-shape sub-windows (see Figure 5) that are determined individually for every pixel in the image, rather than modelling the noise in the overall image. The algorithm is briefly described as follows:
(1)
Select a window of 5 × 5 pixels.
(2)
Calculate the mean n i ( i = 1 , 2 . . . 9 ) and variance α i ( i = 1 , 2 . . . 9 ) of nine masks.
(3)
Select one mask with the lowest variance α k and mean n i , and calculate the weights of every pixel within the k t h mask using the following equation:
ω i = e | Δ i | ,   Δ i = g r a y i n k
(4)
Calculate the output using Equation (28):
g r a y o u t p u t = ( i = 1 M W i g r a y i ) / ( i = 1 M ω i )
where M is the number of pixels in the k t h mask and g r a y i ( i = 1 , 2 . . . 9 ) is the intensity.
With the filter algorithm above, the results of removing the DISP image noise are depicted in Figure 6, which demonstrates the effectiveness of the proposed approach.

3.2. DISP Image Orthorectification and Accuracy Analysis

3.2.1. DISP Image Orthorectification

Since sufficient numbers of GCPs are not observed in each DISP image, TPs are identified to tie overlapping images together and solve for the coefficients of the 2OPE-RM. The study area consists of 355 DISP images (there are 444 DISP images in total, but we only employed 355 high-quality images). Thus, it is impractical to construct a block in the entire study area and then solve for the orthorectification parameters of all DISP images simultaneously because such a large block will produce a significantly large number of observation equations, resulting in a huge computational burden during matrix inversion. Therefore, this paper divides the study area into 24 blocks consisting of various DISP images (see Figure 7a). Each block was rectified independently. For example, Block 1 consists of nine images in Figure 7b, in which 20 GCPs and 29 tie points were identified and measured. The 20 GCPs are employed 36 times to establish 72 observation equations. The 29 TPs are employed 60 times to establish 120 observation equations. Thus, 192 observation equations (72 + 120 = 192) are established in Block 1. There are 166 unknown parameters (9 × 12 + 29 × 2 = 166). There are 26 redundant observations (i.e., 192 − 166 = 26), meeting the requirements of least squares adjustment.
With the 192 observation equations established using Equation (12); the parameters used to rectify the nine DISP images are solved simultaneously using Equation (15). The 2D coordinates of TPs are obtained using Equations (19) and (20). With the solved coefficients and TP coordinates in each image, orthorectification is performed for each DISP image at a GSD of 2.0 m. Figure 8a is part of one orthorectified DISP image.
The computational accuracies of TPs using the 2OPE-RM are evaluated by Equation (23). The standard deviations of TPs ( μ X and μ Y ) are averagely 0.34 m and 0.23 m, respectively. In addition, the “absolute” accuracy of the orthorectified aerial photo created in 1961 is calculated using the following equations:
Δ X R M S E = k = 1 n ( X k x k ) 2 n 1
Δ Y R M S E = k = 1 n ( Y k y k ) 2 n 1
where X k and Y k are XY-coordinates of TPs in the orthorectified DISP image, x k and y k are XY-coordinates in the orthorectified aerial photo created in 1961, and n is the total number of TPs. Using Equations (29) and (30), Δ X R M S E and Δ Y R M S E are 2.0 m and 1.6 m, respectively. These values are equivalent to approximately 2.0 pixels in the orthorectified DISP imagery.

3.2.2. Accuracy Comparison Analysis

Accuracy comparison between the DISP images orthorectified using the traditional second-order polynomial model and the 2OPE-RM was conducted. Two test fields, which are located in mountainous and flat areas, were selected for the accuracy comparison.
(1)
The Bameng Field is a mountainous area located in Bameng County to the west of the city of Baise, Guangxi, China, at 23.671°N to 24.135°N and 106.941°W to 107.698°W. This test area covers the entire DS1106-2119DF107a image. The maximum and minimum elevations are 1128 m and 790 m, respectively, above mean sea level (MSL). Therefore, the relief displacement is significant. There are 12 GCPs and seven TPs scattered throughout the test field. The 12 GCPs are used for second order polynomial equations to solve for the 12 rectification coefficients, and the 12 GCPs and seven TPs are used in the 2OPE-RM to calculate the coefficients. Twenty-three checkpoints were chosen to evaluate the achievable accuracy. The orthorectified aerial photo provided by the Bureau of Guangxi Geomatics and Geographic Information is considered to represent the “true” values for validation. The results are listed in Table 1.
(2)
The Longzhou field is a flat area located in Longzhou County to the west of the city of Chongzuo, Guangxi, China, at 22.105°N to 22.469°N and 106.593°W to 106.878°W. This test field completely covers the entire DS1106-2119DF110a image. In this test field, 11 GCPs and seven TPs are scattered throughout the DISP image. The same GCPs are employed in the traditional second order polynomial model and the 2OPE-RM. Twenty-three checkpoints were chosen to evaluate the accuracy. The planimetric accuracies of the two models relative to the orthorectified aerial image are shown in Table 1.

3.3. Image Mosaicking

Based on the individual image orthorectification above, the next work is to mosaic the individual orthorectified DISP images into an image map. First, the characteristics of the study area and DISP images must be understood:
  • The study area covers 23,790.8 km2 (between latitudes 20.54°N and 26.24°N and longitudes 104.26°E and 112.04°E), which consists of 355 DISP images that total 100 GB. A good mosaicking scheme may save computational time and computer storage;
  • The study area is located in a karst landscape, where mountainous and hilly terrain areas account for two-thirds of the total area;
  • The overlap between neighboring images must be less than 30%; and
  • The study area is covered by five strips of DISP images from four missions (Figure 9).
To minimize the influence of error propagation and avoid repeatedly sampling images, based on the characteristics above, the mosaicking is designed as follows (see Figure 9):
  • The 16 DISP images from Mission 1106 were first mosaicked, covering the western portion of the study area. The mosaicked map is depicted in Figure 10a. Twenty DISP images from Mission 1102-2 were mosaicked, and the mosaicked map is depicted in Figure 10b. Twenty-eight DISP images from Mission 1106 were mosaicked, and the mosaicked map is depicted in Figure 10c. Twenty-three DISP images from Mission 1106 were mosaicked, as the mosaicked map is depicted in Figure 10d. Finally, 18 DISP images from Mission 1106 were mosaicked, and the mosaicked map is depicted in Figure 10d, covering the eastern portion of the area; and
  • With the five mosaicked maps above, a map image of the entire study area was assembled by merging the five mosaicked images. The order of mosaicking is from the east and west to the middle of the study area (see Figure 10).

3.4. Radiance Balance

Due to the differences in the imaging date/time and different imaging conditions during different missions, brightness differences between neighboring strips are unavoidable. In addition, patchwork lines are also unavoidable. To produce a seamless mosaic of the entire study area, this paper used a histogram equalization method to adjust the brightnesses of two neighboring strips. The boundary line was chosen along the center image, and overlapping areas were feathered. Figure 11 shows the result of the radiance balance.

3.5. Mocsaicking Result and Accuracy Evaluation

The entire study area has been mosaicked by 355 orthorectified DISP images (Figure 12a). A mountainous area located in Du’an County (Figure 12b) and a flat area located in Xingbin County (Figure 12c) are select as the samples for accurate validation. Seventy-eight GCPs, which were measured by RTK GPS measurements, are uniformly distributed in other countries throughout the entire study area. These include 25 check points (CPs) scattered throughout the two test fields. Δ X R M S E and Δ Y R M S E in Equations (26) and (27) are used to measure the accuracy. The results are listed in Table 2. As shown in Table 2, the accuracy in flat areas is better than that in mountainous areas, and the overall accuracies of the entire study area are 2.11 m and 1.74 m.

4. Discussions

The results of the accuracy comparison for the DISP images orthorectified by the 2OPE-RM and by traditional second order polynomial model [14] are listed in Table 1. As observed from Table 1, it is demonstrated that the RMSEs of XY-coordinates in the DISP images orthorectified by the 2OPE-RM are smaller than those orthorectified by traditional second order polynomial model in both mountainous and flat areas. With the experimental result, it can be concluded that:
(1)
The proposed 2OPE-RM method can successfully solve the problems below when orthorectifying the DISP images that:
(a)
Each of the original DISP image has insufficient GCPs;
(b)
The camera’s imaging model is unknown; and
(c)
The camera’s interior orientation parameters (IOPs) including camera’s principal point coordinates, focal length, and lens distortion parameters are unknown.
(2)
The proposed 2OPE-RM is capable of obtaining a higher accuracy than the traditional second order polynomial method does when orthorectifying the images under the above conditions (see Table 2).
(3)
Although the proposed 2OPE-RM method is experimented and validated on the DISP images with a satisfied accuracy, it should be suitable for the high-resolution of satellite images whose imaging model and whose camera IOPs are not released.
The major limitations of the proposed 2OPE-RM method include:
(a)
The proposed method needs a lot of tie points, which tie all images together. As observed in Table 1, the accuracy of images rectified by the proposed method can be increased with increasing the number of TPs. For example, if only the 12 GCPs are used in Bameng Field, the RMSEs of the XY-coordinates are 1.96 m and 1.84 m, respectively. If seven TPs are added in addition to 12 GCPs, the RMSEs of XY- coordinates reach 1.85 m and 1.69 m, respectively. The accuracy of the rectification result has been improved. Thus, the more the tie points, the higher the accuracy of orthorectification.
(b)
The proposed method is time-consuming and labor-intensive, because a lot of tie points, which are usually feature points in images, are manually selected and measured. Although a semi-automation of measurement and selection of TPs are used in this paper, a zoom-in window operation for high-accuracy of location of the TPs is usually employed.

5. Conclusions

This paper presents a highly effective, simple, practical mathematical model for the orthorectification of CORONA DISP images from the 1960s, whose interior and exterior parameters are unknown and in which GCPs are lacking. The model is called the second order polynomial equation-based block rectification model (2OPE-RM). With the proposed model, all images can be orthorectified at an accuracy level of 2.0 pixels, corresponding to approximately 2.0–4.0 m with respect to the WGS 84 datum. All of the images covering the entire karst area of Guangxi, China, are assembled into a high-quality image map. The sampled distance of the assembled mosaicking map is 2.0 m. The proposed model can solve the problems associated with the traditional second order polynomial model, such as lack of GCPs, yielding acceptable and improved accuracy. The assembled image map of the entire rock desertification area in Guangxi, China, will be delivered to the Guangxi Geological Library for use by the research community.

Acknowledgments

This paper is financially supported by the National Natural Science of China under Grant numbers 41431179 and 41162011, Guangxi Natural Science Foundation under grant numbers 2015GXNSFDA139032, Guangxi Science & Technology Development Program under the contract number GuiKeHe 14123001-4, and GuangXi Key Laboratory of Spatial Information and Geomatics Program (Contract No. GuiKeNeng 140452401, 140452409, 151400701 and 151400712), the “BaGuiScholars” program of the provincial government of Guangxi.

Author Contributions

Guoqing Zhou contributes most to this manuscript. He has made the main contribution to the whole idea and writing the initial draft for this manuscript. Tao Yue contributed to re-organization of the framework and the multiple reviewing of the manuscript. Yujun Shi supported data processing and analysis. Rongting Zhang and Jingjin Huang contributed to editing and reviewing the manuscript.

Conflicts of Interest

The authors declare no conflict of interest.

Abbreviations

The following abbreviations are used in this manuscript:
GCP
Ground Control Point
DISP
Declassified Intelligence Satellite Photography
KRD
Karst Rocky Desertification
IOP
Interior Orientation Parameter
EOP
Exterior Orientation Parameter
RPC
Rational Polynomial Coefficient
2OPE-RM
Second Order Polynomial Equation-Based Rectification Model
GSD
Ground Sampled Distance
TP
Tie Point
CP
Check Point
RTK
Real–Time Kinematic
SPOT
Systeme Probatoire d’Observation de la Terre
SAR
Synthetic Aperture Radar

References

  1. Zhou, G.; Huang, J.; Tao, X.; Luo, Q.; Zhang, R.; Liu, Z. Overview of 30 years of research on solubility trapping in Chinese karst. Earth-Sci. Rev. 2015, 146, 183–194. [Google Scholar] [CrossRef]
  2. Xu, S.; Jiang, Z. Preliminary estimate the source and sink relation of karstification and atmospheric greenhouse gas CO2. Chin. Sci. Bull. 1997, 42, 953–956. [Google Scholar]
  3. Jiang, Z.C.; Qin, X.Q. Calculation of atmospheric CO2 sink formed in karst progresses of the karst divided regions in China. Geosci. Sin. 2011, 30, 363–367. [Google Scholar]
  4. Zhou, G.; Jezek, K.C. Satellite photograph mosaics of Greenland from the 1960s era. Int. J. Remote Sens. 2002, 23, 1143–1159. [Google Scholar] [CrossRef]
  5. Zhou, G.; Jezek, K.C.; Wright, W.; Granger, J. Orthorectification of 1960s satellite photographs covering Greenland. IEEE Trans. Geosci. Remote Sens. 2002, 40, 1247–1259. [Google Scholar] [CrossRef]
  6. Kim, K.; Jezek, K.C.; Sohn, H. Ice shelf advance and retreat rates along the coast of Queen Maud Land, Antarctica. J. Geophys. Res. 2001, 10, 7097–7106. [Google Scholar] [CrossRef]
  7. Kim, K.; Jezek, K.C.; Liu, H. Orthorectified image mosaic of Antarctica from 1963 Argon satellite photography: Image processing and glaciological applications. Int. J. Remote Sens. 2007, 28, 5357–5373. [Google Scholar] [CrossRef]
  8. Tao, C.V.; Hu, Y. A comprehensive study of the rational function model for photogrammetric processing. Photogramm. Eng. Remote Sens. 2001, 67, 1347–1357. [Google Scholar]
  9. Tao, C.V.; Hu, Y. 3D reconstruction methods on the rational function model. Photogramm. Eng. Remote Sens. 2002, 68, 705–714. [Google Scholar]
  10. Yang, X.H. Accuracy of rational function approximation in photogrammetry. In Proceedings of the ASPRS Annual Conference, Amsterdam, The Netherlands, 16–22 July 2000; pp. 22–26.
  11. Liu, J.; W, D.H.; Mao, G.M. High accuracy stereo positioning of IKONOS satellite image based on RPC model. Bull. Surv. Mapp. 2004, 0911, 1–4. [Google Scholar]
  12. Huang, G.M.; Yue, X.J.; Zhao, Z.; Fan, H.D. Block adjustment with airborne SAR images based on polynomial ortho-rectification. Geomat. Inf. Sci. Wuhan Univ. 2008, 33, 569–572. [Google Scholar]
  13. Grodecki, J.; Dial, G. Block adjustment of high-resolution satellite images described by rational polynomials. Photogramm. Eng. Remote Sens. 2003, 69, 59–68. [Google Scholar] [CrossRef]
  14. Zoej, M.J.V.; Mansourian, A.; Mojaradi, B.; Sadeghian, S. 2D geometric correction of IKONOS imagery using genetic algorithm. In Proceedings of the ISPRS Commission IV, Symposium 2002 Geospatial Theory, Processing and Applications, Ottawa, ON, Canada, 9–12 July 2002; pp. 9–12.
Figure 1. The imaging geometry of declassified intelligence satellite photography (DISP) from the CORONA mission.
Figure 1. The imaging geometry of declassified intelligence satellite photography (DISP) from the CORONA mission.
Remotesensing 08 00680 g001
Figure 2. Illustration of the second-order polynomial equation-based rectification model (2OPE-RM).
Figure 2. Illustration of the second-order polynomial equation-based rectification model (2OPE-RM).
Remotesensing 08 00680 g002
Figure 3. Study area, which is located in Guangxi, China, with the encompassing 23,790.8 km2.
Figure 3. Study area, which is located in Guangxi, China, with the encompassing 23,790.8 km2.
Remotesensing 08 00680 g003
Figure 4. DISP image dataset (There are 444 DISP images from five orbits of different missions covering the whole study area).
Figure 4. DISP image dataset (There are 444 DISP images from five orbits of different missions covering the whole study area).
Remotesensing 08 00680 g004
Figure 5. The adaptive filter algorithm.
Figure 5. The adaptive filter algorithm.
Remotesensing 08 00680 g005
Figure 6. The results of noise removing: (a) The original image, and (b) The filtered image.
Figure 6. The results of noise removing: (a) The original image, and (b) The filtered image.
Remotesensing 08 00680 g006
Figure 7. The polynomial block adjustment: (a) Design of polynomial block adjustment in the entire area, which is divided into 24 blocks; (b) One block, which is used to explain establishment of the observation equations.
Figure 7. The polynomial block adjustment: (a) Design of polynomial block adjustment in the entire area, which is divided into 24 blocks; (b) One block, which is used to explain establishment of the observation equations.
Remotesensing 08 00680 g007
Figure 8. Accuracy verification using orthorectified aerial photos from 1961: (a) DISP image orthorectified using the proposed method; and (b) orthorectified aerial photo from 1961.
Figure 8. Accuracy verification using orthorectified aerial photos from 1961: (a) DISP image orthorectified using the proposed method; and (b) orthorectified aerial photo from 1961.
Remotesensing 08 00680 g008
Figure 9. Mosaicking process from east/west to the middle.
Figure 9. Mosaicking process from east/west to the middle.
Remotesensing 08 00680 g009
Figure 10. The mosaicked images of the various missions.
Figure 10. The mosaicked images of the various missions.
Remotesensing 08 00680 g010
Figure 11. Radiometric balance of neighboring strips: (a) The image before radiometric balancing; (b) The image after radiometric balancing.
Figure 11. Radiometric balance of neighboring strips: (a) The image before radiometric balancing; (b) The image after radiometric balancing.
Remotesensing 08 00680 g011
Figure 12. Mosaic results: (a) a completely assembled 1960s mosaic of the Guangxi karst area; (b) accuracy validation in a mountainous area; and (c) accuracy validation in a flat area.
Figure 12. Mosaic results: (a) a completely assembled 1960s mosaic of the Guangxi karst area; (b) accuracy validation in a mountainous area; and (c) accuracy validation in a flat area.
Remotesensing 08 00680 g012
Table 1. Accuracy comparison of mountainous and flat areas.
Table 1. Accuracy comparison of mountainous and flat areas.
Test AreaBameng Field
(Located in A Mountainous Area)
Longzhou Field
(Located in A Flat Area)
Models Δ X R M S E (m) Δ Y R M S E (m) Δ X R M S E (m) Δ Y R M S E (m)
Second order polynomial model1.961.841.851.57
Our model1.851.691.671.49
Table 2. Final accuracies of the assembled DISP image map in the study area.
Table 2. Final accuracies of the assembled DISP image map in the study area.
Area Δ X R M S E (m) Δ Y R M S E (m)
Mountainous area2.071.60
Flat area1.861.79
The entire study area2.111.74

Share and Cite

MDPI and ACS Style

Zhou, G.; Yue, T.; Shi, Y.; Zhang, R.; Huang, J. Second-Order Polynomial Equation-Based Block Adjustment for Orthorectification of DISP Imagery. Remote Sens. 2016, 8, 680. https://0-doi-org.brum.beds.ac.uk/10.3390/rs8080680

AMA Style

Zhou G, Yue T, Shi Y, Zhang R, Huang J. Second-Order Polynomial Equation-Based Block Adjustment for Orthorectification of DISP Imagery. Remote Sensing. 2016; 8(8):680. https://0-doi-org.brum.beds.ac.uk/10.3390/rs8080680

Chicago/Turabian Style

Zhou, Guoqing, Tao Yue, Yujun Shi, Rongting Zhang, and Jingjin Huang. 2016. "Second-Order Polynomial Equation-Based Block Adjustment for Orthorectification of DISP Imagery" Remote Sensing 8, no. 8: 680. https://0-doi-org.brum.beds.ac.uk/10.3390/rs8080680

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop