Assumptions:
The height above sea level in metres is derived from an elevation image. The elevation for each GCP can be derived from an elevation channel, corrected to UTM.
NGCP - number of GCPs
P(i) - pixel number in the input SAR image for the (i)th GCP
L(i) - line number in the input SAR image for the (i)th GCP
E(i) - easting of the pixel number for the (i)th GCP
N(i) - northing of the line number for the (i)th GCP
H(i) - height (elevation) above sea level for the (i)th GCP
Initial estimates of the flight line parameters are given as
follows:
ALTI - height of the aircraft above sea level HEAD - heading of the aircraft, clockwise from due north E0 - easting of a point flown over by the aircraft (POINT(1)) N0 - northing of a point flown over by the aircraft (POINT(2))Estimating Ground Range Distance for each GCP:
FLIGHT calculates two estimates of the ground range distance of each GCP from the flight path. If the flight path parameters are accurate, the values for the two estimates should be the same. FLIGHT adjusts the values of the flight path parameters until the sum of the squares of the differences between the two estimates is less than some threshold value (MAPTOL).
First Estimate for Ground Range:
The first method of estimating ground range uses the geometry of the flight path, and no elevation data. The flight path is projected onto a flat earth surface. The ground range is the length of the perpendicular line connecting the GCP to the flight path. Given the easting and northing of one point (P) on the flight path, the easting and northing of the GCP, and a heading (which is the angle between the flight path and due north, and also the angle (theta) between the perpendicular line connecting the GCP to the flight path and due west), it is possible to calculate the coordinates of the intercept point (P2) between the flight path and the perpendicular line to the GCP. The ground range is the absolute distance between the intercept point on the flight path and the GCP. See the diagram below.
N
^
| ^ - Flight Path
| /
| /
| x P(E0,N0)
| /
| /
| o<----P2 ( EI(i),NI(i) )
| / \
Heading-->|--/ \<--Ground Range
Angle | / \ (perpendicular to
|/theta>|\ flight path)
W<------------------------+---------o---------------
/| |
/ | GCP ( E(i),N(i) )
/ |
/ |
/ |
/ |
/ |
/ |
^ |
N = northing; E = easting
The equation for the flight path is as follows:
N = N0 - tan(90-HEAD)*E0 + tan(90-HEAD)*E
The equation for the perpendicular line between the GCP and the
flight path is as follows:
N = N(i) + tan(HEAD)*E(i) - tan(HEAD)*E
Solving for the interception point for these two lines gives the
following:
EI(i) = easting of the intercept point for the (i)th GCP
NI(i) = northing of the intercept point for the (i)th GCP
N(i) - N0 + tan(HEAD)*E(i) + tan(90-HEAD)*E0
EI(i) = --------------------------------------------
tan(HEAD) + tan(90-HEAD)
NI(i) = N0 + tan(90-HEAD)*(EI(i) - E0)
The above equations cannot be used if HEAD is a multiple of 90.If HEAD = 0 (north) or 180 (south), then the intercept point is as follows:
EI(i) = E0
NI(i) = N(i)
If HEAD = 90 (east) or 270 (west), then the intercept point is as
follows:
EI(i) = E(i)
NI(i) = N0
The first ground range distance estimate for GCP(i) is:
G1(i) = square root of ((EI(i) - E(i))**2 + (NI(i) - N(i))**2)
Second Estimate for Ground Range:The second method of estimating ground range uses the altitude of the aircraft above the GCP and the slant range distance to the GCP. If the input SAR data is a ground range image, then the ground range distance to the GCP is converted to slant range distance, using the assumed aircraft height above the ground. This conversion is done before estimating the true ground range distance, given the altitude of the aircraft and GCP above sea level.
DELAY = time delay to first input SAR pixel
C = speed of light (299.793 metres per microsecond)
S0 = slant range distance to first input SAR pixel
G0 = ground range distance to first input SAR pixel
INPXSZ(1) = input SAR pixel spacing in range direction
HEIGHT = assumed aircraft height above ground (if input
SAR image is ground range data)
S0 = DELAY * C / 2
G0 = SQRT(S0**2 - HEIGHT**2)
If input SAR data is a slant range image (RANGETYP="SLANT"):Calculate the slant range distance, given pixel number P(i):
S(i) = S0 + (P(i)-1)*INPXSZ(1)
The second ground range distance estimate for GCP(i) is:
G2(i) = SQRT(S(i)**2 - (ALTI - H(i))**2)
If input SAR data is a ground range image (RANGETYP="GROUND"):Calculate the ground range distance, given pixel number P(i):
G(i) = G0 + (P(i)-1)*INPXSZ(1)
The second ground range distance estimate, taking elevation data
into account for for GCP(i) is:
G2(i) = SQRT(G(i)**2 + HEIGHT**2 - (ALTI - H(i))**2)
Then, the estimation error is calculated as follows:
SUM for i = 1 to n of (( G1(i) - G2(i) )**2 )/n
ERROR**2 = sum of ((G1(i) - G2(i))**2) over all GCPs
divided by number of GCPs
If ERROR is less than or equal to MAPTOL (error threshold), then the
flight path parameters are considered accurate, and the FLIGHT
program terminates and saves these parameters in the parameter file.Adjusting Flight Path Parameters:
If ERROR is greater than MAPTOL, then the flight path parameters are adjusted, and the error is calculated again until the error threshold is met. The value of ERROR is calculated for six sets of adjusted flight path parameters. The following list contains adjustment offsets to the three flight path parameters:
DALTI = adjustment offset for altitude
DHEAD = adjustment offset for heading
DPOINT = adjustment offset for point on flight path
There are six different adjustments which are attempted. Of the six
adjustments, the one with the minimum ERROR is chosen. If ERROR is
less than MAPTOL, then FLIGHT outputs the flight parameters to the
parameter file and terminates. Otherwise, another iteration of six
cases is calculated, using the best previous case as the base.Given the above adjustment offsets to the flight path parameters, the six adjustments are as follows (flight parameters which are not mentioned remain unchanged):
Case 1. Add DALTI to altitude.
Case 2. Subtract DALTI from altitude.
Case 3. Add DHEAD to heading.
Case 4. Subtract DHEAD from heading.
Case 5. If heading is between 0-45, 135-225 or 315-360, then
add DPOINT to easting of point on flight path.
If heading is anything else, then
add DPOINT to northing of point on flight path.
Case 6. If heading is between 0-45, 135-225 or 315-360, then
subtract DPOINT from easting of point on flight path.
If heading is anything else, then
subtract DPOINT from northing of point on flight path.
The values for DALTI, DHEAD and DPOINT are derived from MAPTOL.
DALTI = MAPTOL
DHEAD = arctan(MAPTOL / length of input SAR image)
where the length of the input SAR image equals the number of lines
in the uncorrected database, times the azimuth pixel spacing
(INPXSZ(2)).
DPOINT = MAPTOL
Normally, the value of ERROR decreases with each iteration. However,
if the ERROR for the current iteration is not less than the value of
ERROR for the previous iteration, then DALTI, DHEAD, and DPOINT are
divided in half for the next iteration. If the value of ERROR is
less than the previous value of ERROR, then the original values for
DALTI, DHEAD, and DPOINT are restored for the next iteration. This
method handles the situation where the algorithm does not converge.
Iterations are repeated until the algorithm converges and ERROR is
less than MAPTOL, or until a total of 500 iterations have been
performed. Then, the program terminates, and the report is produced.Determining Coefficients of Line Polynomial Function:
Once values for OALTI, OHEAD, and OPOINT are determined accurately, then the coefficients of the polynomial equation which calculates the input SAR line numbers, given their distance along the flight path from the OPOINT coordinate, can be calculated. Coefficients are determined using least squares regression analysis.
The intercept point of a perpendicular line connecting each GCP to flight path has been previously calculated. The distance between OPOINT and the intercept point for the (i)th GCP is calculated as follows:
D(i) = distance from OPOINT to intercept point for (i)th GCP
EI(i) = easting of intercept point for (i)th GCP
NI(i) = northing of intercept point for (i)th GCP
E0 = easting of point of flight path (OPOINT(1))
N0 = northing of point of flight path (OPOINT(2))
D(i) = SQRT((EI(i) - E0)**2 + (NI(i) - N0)**2)
The sign of D(i) is positive if NI(i) > N0.The sign of D(i) is negative if NI(i) < N0.
Given values for distance (D(i)) and line number (L(i)) for each GCP, regression analysis is used to determine a polynomial equation of user specified ORDER which best fits these points.
L = a + (b * D ) + (c * D**2) + (d * D**3) + (e * D**4)
+ (f * D**5) + (g * D**6) + (h * D**7) + (i * D**8)