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

Guoqing Zhou 1,*, Tao Yue 1,*, Yujun Shi 2, Rongting Zhang 3 and Jingjin Huang 3 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; shiyujun1988@163.com 3 The Center for Remote Sensing, Tianjin University, No. 92, Weijin Road, Nankai District, Tianjin 300072, China; zrt65@tju.edu.cn (R.Z.); jingjin_huang@163.com (J.H.) * Correspondence: gzhou@glut.edu.cn (G.Z.); yuetao@glut.edu.cn (T.Y.); Tel.: +86-773-589-6073 (G.Z.); +86-773-589-1963 (T.Y.)


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 CO 2 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 km 2 , accounting for 19.8% of the total KRD area in China in 2005.It shrunk to 19,260.00km 2 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.

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 − WVU(W V U ) is a camera coordinate system, o − xy is an image plane system, and O − XYZ 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.
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: 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.

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]:  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: where x and y are image coordinates; ∆x and ∆y are image distortions in the xand 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.

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]: where α = (a 0 , a 1 , a 2 , a 3 , a For a given GCP, Equation (3) can be linearized using a Taylor series and is expressed as follows: 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).
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: where: The least-squares estimation, i.e., V T PV = min, gives the solutions of the coefficients of the second order polynomial equation below: We can further obtain the following expressions from Equation (7): where a 0 i , b 0 i are initial values; ∆a j i , ∆b j i are increases during each iteration; and N ite is the number of iterations.

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: Then, Equation ( 10) can be rewritten as follows: where: 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: where: The 1 st image with N 1 GCPs The N th image with N n GCPs and: The 1 st image with M 1 TPs The i th image with M i 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: With least-squares estimation, the normal equation matrix can be written as follows: Thus, the solution of the unknown parameters is given by Equation ( 15): where Q ij (i, j = 1, 2) gives the components of the covariance matrix, which is the inverse of the normal matrix, as shown in Equation ( 16): The coefficients of the 2OPE-RM in each image and the 2D coordinates (XY) of each TP are as follows: where X i , Y i are coordinates of the i-th TP in image t i ; ∆X i , ∆Y i are increases in X i and Y i ; a 0 i and b 0 i 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): 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: To evaluate the accuracies of TPs, assuming that there are n TPs, the average of δ X i is as follows: where n is the number of TPs.
Remote Sens. 2016, 8, x FOR PEER 7 of 18 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): 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: To evaluate the accuracies of TPs, assuming that there are n TPs, the average of i X  is as follows: where n is the number of TPs.

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:

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:

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 ul , Y ul ) , (X ll , Y ll ) , (X ur , Y ur ) , (X lr , Y lr ) .
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 ul , X ll , X ur , X lr ) , X max = min(X ul , X ll , X ur , X lr ) • Determination of the resulting image's size: The size of the resulting image can be determined by M and N as follows: where M = row, N = col, and Y GSD , X GSD 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 ori , y ori ) and ground coordinates X gro , Y gro , the ground coordinates should be transformed into the coordinates of the resulting image (x re , y re ) as follows: where Y gro , X gro 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 ori via bilinear resampling interpolation.

•
Assign the gray value g ori to point P as g res 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.

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.8km 2 .The study area is in the south central subtropics of China.Aerial photos: Five aerial photos with film formats of 18 18  cm 2 from 1961 were acquired at a photographic scale of 1:14,000.Each photo covers approximately 6.35 km 2 .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.• 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 cm 2 from 1961 were acquired at a photographic scale of 1:14,000.Each photo covers approximately 6.35 km 2 .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. 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  cm 2 from 1961 were acquired at a photographic scale of 1:14,000.Each photo covers approximately 6.35 km 2 .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.

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.
(3) Select one mask with the lowest variance α k and mean n i , and calculate the weights of every pixel within the k th mask using the following equation: (4) Calculate the output using Equation (28): where M is the number of pixels in the k th mask and gray 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.

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 of nine masks.
(3) Select one mask with the lowest variance k  and mean i n , and calculate the weights of every pixel within the th k mask using the following equation: (4) Calculate the output using Equation (28): where M is the number of pixels in the th k mask and 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.

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 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:

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 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: 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 RMSE and ∆Y RMSE are 2.0 m and 1.6 m, respectively.These values are equivalent to approximately 2.0 pixels in the orthorectified DISP imagery.
where k X and k Y are XY-coordinates of TPs in the orthorectified DISP image, k     )

Accuracy Comparison Analysis
where

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

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: 1.
The study area covers 23,790.8km 2 (between latitudes 20.54The study area is located in a karst landscape, where mountainous and hilly terrain areas account for two-thirds of the total area; 3.
The overlap between neighboring images must be less than 30%; and 4.
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): 1.
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 2.
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).
Remote Sens. 2016, 8, x FOR PEER 14 of 18 2. 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).

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.2. 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).

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.

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.

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.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.

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 RMSE and ∆Y RMSE 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.

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:

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: (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.

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.

Figure 1 .
Figure 1.The imaging geometry of declassified intelligence satellite photography (DISP) from the CORONA mission.
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.

Figure 1 .
Figure 1.The imaging geometry of declassified intelligence satellite photography (DISP) from the CORONA mission.

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 xand y-directions, respectively; and X and Y are 2D coordinates in a given map coordinate system.

Figure 3 .
Figure 3. Study area, which is located in Guangxi, China, with the encompassing 23,790.8km 2 .

Figure 4 .
Figure 4. DISP image dataset (There are 444 DISP images from five orbits of different missions covering the whole study area).

Figure 3 .
Figure 3. Study area, which is located in Guangxi, China, with the encompassing 23,790.8km 2 .

Figure 3 .
Figure 3. Study area, which is located in Guangxi, China, with the encompassing 23,790.8km 2 .

Figure 4 .
Figure 4. DISP image dataset (There are 444 DISP images from five orbits of different missions covering the whole study area).

Figure 4 .
Figure 4. DISP image dataset (There are 444 DISP images from five orbits of different missions covering the whole study area).

Figure 5 .
Figure 5.The adaptive filter algorithm.Figure 5.The adaptive filter algorithm.

Figure 5 .
Figure 5.The adaptive filter algorithm.Figure 5.The adaptive filter algorithm.

Figure 6 .
Figure 6.The results of noise removing: (a) The original image, and (b) The filtered image.

Figure 6 .
Figure 6.The results of noise removing: (a) The original image, and (b) The filtered image.
x and k y are XY-coordinates in the orthorectified aerial photo created in 1961, and n is the total number of TPs.Using Equation (29) and Equation (30), and 1.6 m, respectively.These values are equivalent to approximately 2.0 pixels in the orthorectified DISP imagery.

Figure 7 .
Figure 7.The polynomial block adjustment: (a) Design of polynomial block adjustment in the area, which is divided into 24 blocks.(b) One block, which is used to explain establishment of the observation equations.

Figure 8 .
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 7 .
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.

Y
are XY-coordinates of TPs in the orthorectified DISP image, in the orthorectified aerial photo created in 1961, and n is the total number of TPs.Using Equation (29) and Equation (30), and 1.6 m, respectively.These values are equivalent to approximately 2.0 pixels in the orthorectified DISP imagery.

Figure 7 .
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 8 .
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 .
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 9 .
Figure 9. Mosaicking process from east/west to the middle.

Figure 10 .
Figure 10.The mosaicked images of the various missions.

Figure 9 .
Figure 9. Mosaicking process from east/west to the middle.

Figure 9 .
Figure 9. Mosaicking process from east/west to the middle.

Figure 10 .
Figure 10.The mosaicked images of the various missions.

Figure 10 .
Figure 10.The mosaicked images of the various missions.

Figure 11 .
Figure 11.Radiometric balance of neighboring strips: (a) The image before radiometric balancing, (b) The image after radiometric balancing.
26) and (27) are used to measure the accuracy.The results are listed in Table

Figure 11 .
Figure 11.Radiometric balance of neighboring strips: (a) The image before radiometric balancing; (b) The image after radiometric balancing.

Figure 12 .
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.
(a) Each of the original DISP image has insufficient GCPs; (b) The camera's imaging model is unknown; and

Figure 12 .
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.

( 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.

Table 1 .
• 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 Table1.(2)TheLongzhou 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.Accuracy comparison of mountainous and flat areas.
• 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; 2.

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.