422 |
H.-E. Andersen |
m32 = − sin(ω) cos(ϕ) m33 = cos(ω) cos(ϕ)
and (Xo1 , Yo1 , Zo1 ) is the exposure station coordinate corresponding to image 1 (see Fig. 10.2).
In stereo imagery (i.e. when the same feature is imaged in two photographs taken from different perspectives), then two collinearities are formed for each scene point and these rays can be intersected to determine the 3D position of the scene point. Thus, when a feature is imaged within the overlap area on two stereo images, the pair of Eqs. (10.1) and (10.2) can be formed for each image, giving a system of four equations and three unknowns (Xp , Yp , Zp ). This coordinate can therefore be determined via a least-squares solution. As an alternative to the least-squares solution and for the purposes of demonstration, the coordinates of an object on the ground can be obtained using space intersection via the following approach described in [56]. If the image point in each photo (xp , yp ) is described in terms of the coordinate system of the tilted photography (x, y, −f ) and another untilted image coordinate system whose respective axes (x, y, z) are parallel to the axis of the ground coordinate system (X, Y, Z), then the coordinates of object point (Xp , Yp , Zp ) can be expressed as a scaled version of the untilted image 1 coordinate system (xp1 , yp1 , zp1 ):
Xp = λp1 xp1 + Xo1 |
|
|
Yp = λp1 yp |
1 + Yo1 |
(10.4) |
Zp = λp1 zp |
1 + Zo1 |
|
and the untilted coordinates of the image point p can be obtained from the tilted image coordinates via:
xp1 = m111 xp1 + m211 yp1 + m311 zp1 |
|
|
yp |
1 = m121 xp1 + m221 yp1 + m321 zp1 |
(10.5) |
zp |
1 = m131 xp1 + m231 yp1 + m331 zp1 |
|
Since the object point with coordinates (Xp , Yp , Zp ) is imaged in both photos, the coordinate can also be expressed as a scaled version of the untilted image 2 coordinate system (xp2 , yp2 , zp2 ) using Eqs. (10.4).
Solving for the scaling factor λp1 in these six equations yields:
y |
|
(XO |
1 − |
XO |
) |
− |
x |
(YO |
1 − |
YO |
) |
|
|||
λp1 = |
|
p1 |
|
|
|
2 |
|
p1 |
2 |
|
(10.6) |
||||
|
|
|
xp |
1 |
yp |
|
− |
xp |
yp |
|
|
|
|||
|
|
|
|
|
|
1 |
2 |
2 |
|
|
|
|
|||
which can then be used in Eqs. (10.4) to calculate the coordinate of the point (Xp , Yp , Zp ).
Example Given two images taken with a digital camera with a pixel size of 5.3 microns and a calibrated focal length of 35.1138 mm, and with exterior orientation parameters (X, Y, Z, ω, ϕ, κ ) of (404829.1, 7029981.0, 1202.1, 5.71, 2.42, 93.14) for image 1 and (404962.9, 7029992.0, 1197.9, 4.65, 4.87, 91.63) for image 2, calcu-
10 High-Resolution Three-Dimensional Remote Sensing for Forest Measurement |
423 |
late the object space coordinate for a feature with pixel coordinates (3445.88, 2435.13) in image 1 and (3561.62, 1482.28) in image 2.
Solution Using Eqs. (10.3), the elements of the rotation matrix for each image can be calculated from the (ω, ϕ, κ ) angles provided in the exterior orientation. This yields:
|
|
|
m11 = |
−0.055(image1) |
|
|
−0.029(image2) |
|
|
|
|
||||||||||||||
|
|
|
m12 = |
0.99(image1) |
|
|
1.0(image2) |
|
|
|
|
|
|||||||||||||
|
|
|
m13 = |
0.10(image1) |
|
|
0.083(image2) |
|
|
|
|
|
|||||||||||||
|
|
|
m21 = |
−0.98(image1) |
|
|
−1.0(image2) |
|
|
|
|
|
|||||||||||||
|
|
|
m22 = |
−0.059(image1) |
|
|
−0.035(image2) |
|
|
|
|
||||||||||||||
|
|
|
m23 |
= |
−0.047(image1) |
|
|
−0.087(image2) |
|
|
|
|
|||||||||||||
|
|
|
m31 |
= |
0.042(image1) |
|
|
0.085(image2) |
|
|
|
|
|
||||||||||||
|
|
|
m32 |
= |
−0.099(image1) |
|
|
−0.08(image2) |
|
|
|
|
|
||||||||||||
|
= − |
m33 |
= |
0.99(image1) |
|
|
0.99(image2) |
|
= |
|
|
||||||||||||||
p1 |
0.055)(6.90) |
+ − |
|
− |
5.36) |
+ |
(0.042)( |
− |
35.11) |
3.49 |
|||||||||||||||
x |
= |
( |
|
|
( 1.0)( |
|
|
|
|
|
|
|
|||||||||||||
p1 |
(0.99)(6.90) |
+ − |
0.059)( |
− |
5.36) |
+ |
( |
− |
|
− |
|
|
= |
10.66 (10.7) |
|||||||||||
y |
= |
|
( |
|
− |
+ |
|
0.099)( |
35.11) |
|
|
||||||||||||||
p1 |
(0.10)(6.90) |
+ − |
0.047)( |
5.36) |
|
|
|
− |
35.11) |
= − |
33.96 |
||||||||||||||
z |
|
|
( |
|
|
|
|
|
(0.99)( |
|
|
|
|||||||||||||
The untilted image coordinates for image 2 can be calculated in a similar manner, yielding xp2 = −2.88, yp2 = 10.33, zp2 = −34.22. The scaling factor λp1 can then be calculated as:
|
y |
|
(XO |
1 − |
XO ) |
− |
x |
(YO |
− |
YO ) |
|||
λp1 = |
|
p1 |
|
|
|
2 |
p1 |
1 |
2 |
||||
|
|
|
xp |
1 |
yp |
− |
xp |
yp |
2 |
|
|
||
|
|
|
|
|
|
1 |
2 |
|
|
|
|||
= 10.66(404829.1 − 404962.9) − 3.49(7029981.1 − 7029991.7) (3.49)(10.66) − (−2.88)(10.33)
= 21.15
The object coordinates for the feature of interest can then be calculated as:
Xp = λp1 xp1 + Xo1 = (21.15)(3.49) + 404829.1 = 404975.1 Yp = λp1 yp1 + Yo1 = (21.15)(10.66) + 7029981.1 = 7029868.8 Zp = λp1 zp1 + Zo1 = (21.15)(−35.11) + 1202.1 = 484.1
With the advent of the digital computer, it is possible to measure tree heights on digital imagery within a digital (or softcopy) photogrammetric workstation environment. In a digital photogrammetric workstation, a pair of overlapping photographs can be viewed in stereo on a computer monitor using either special hardware (shutter glasses and an emitter to synchronize the display and the glasses) or color anaglyph display and glasses with red-cyan filters (see Fig. 10.3).
Furthermore, a digital photogrammetric workstation provides the capability to digitize features within a stereo model, such as the 3D coordinate of a tree top and
424 |
H.-E. Andersen |
Fig. 10.3 Stereo aerial photography of a forest scene viewed in a digital photogrammetric environment using color anaglyph, upper Tanana valley of interior Alaska, USA. (Red-blue glasses are necessary to view this scene in stereo). Significant radial displacement: apparent shift of an object having height in relation to its base in an image with a central projection of trees (layover) is evident in the top left corner of the scene
crown base (tree height) or a (planimetrically-correct) polygon delineating a distinct forest condition class, such as its size class2 or species class3 or density class.4
As the collinearity conditions indicate above, it is essential to know the coordinates of each camera station and the elements of the rotation matrix (M) of each camera before aerial photographs can be used to acquire three-dimensional measurements of forest features (Xp , Yp , Zp ). Typically, this information is obtained through the use of ground control points, which are features on the ground, with known (Xp , Yp , Zp ) coordinates, that are also visible in the overlap area of the stereo pair of aerial photographs (technically, three vertical control points for leveling the model and two horizontal control points for scaling the model is the minimum requirement for controlling a single stereo pair [56]). If additional ground control points are available, a least-squares solution can be obtained with an estimate of uncertainty. When there are a large number of overlapping stereo models covering an area, the requirement of three control points per stereo model is relaxed, and the exterior orientation parameters for all photos within the block can be obtained through a procedure known as a bundle block adjustment [56]. As stated
2Size class refers to the predominant size (diameter, or DBH) of the trees within a stand; e.g. regeneration (DBH < 12.7 cm), poletimber (12.7 cm < DBH < 30.48 cm), sawtimber (DBH > 30.48 cm).
3Species class refers to the predominant species (or species mixture) of the stand (e.g. black spruce, mixed spruce-hemlock, etc.).
4Density class refers to the number of tree stems in a given area (i.e. trees per hectare).
10 High-Resolution Three-Dimensional Remote Sensing for Forest Measurement |
425 |
in [56], a bundle adjustment is a procedure for simultaneously “adjusting all photogrammetric measurements to ground control values in a single solution.” However, when using aerial photographs in a sampling mode, where plots are widely spaced, the need for sufficient ground control can introduce a significant and potentially prohibitive additional cost to the aerial photo-based inventory.
If smaller scale (i.e. lower resolution), controlled aerial photography is available that covers the same area as the uncontrolled large scale photography, ground control points can be measured photogrammetrically in the small scale photography and subsequently used to control the large scale photography (a method known as bridging control) [46]. However, in very remote areas, such as interior Alaska, it is unlikely that even recent small-scale, controlled imagery will be available. Fortunately, in recent years, technology has become available that allows for precise, and accurate, measurement of the position and orientation of the camera in the aircraft at the moment of exposure, which significantly reduces, or even eliminates entirely, the need for surveyed ground control. The use of two tightly-coupled technologies: (1) the global positioning systems (GPS) and (2) inertial measurement unit (IMU), now allow the exterior orientation parameters for each camera station to be obtained without the need for ground control, an approach known as direct georeferencing [34]. The GPS instrument uses a system of satellites to triangulate the position of the camera, while the inertial measurement unit uses a system of accelerometers and gyroscopes to determine the orientation of the camera. Furthermore, because the GPS acquires accurate, but relatively noisy, positional information, while the IMU provides trajectory and orientation information with relatively little noise but with systematic drift, the positional error can be dramatically reduced by merging these two complementary sources of positional information via a Kalman filter signal processing procedure. In the Kalman filter, the estimate of the position at time k + 1 is given by the so-called state estimate equation:
xk+1 = Axk + Buk + wk
yk = Cxk + zk
Kk = APk CT CPk CT + Sz −1
xˆ k+1 = (Axˆ k + Buk ) + Kk (yk+1 − Cxˆ k )
Pk+1 = APk AT + Sw − APk CT S−z 1CPk AT
In the above equation, A, B, and C are matrices that describe how the state changes and can be measured, k is the time index, x is the state of the system, u is the known input to the system, y is the measured output, w is the process noise, z is the measurement noise, Kk is the Kalman gain, Sw is the process noise covariance: Sw = E(wk wTk ), and Sz is the measurement noise covariance: Sz = E(zk zTk ), and P is the estimation error covariance. The first term in the fourth equation is basically A times the estimated position xˆ at time k, plus B times the known input u (IMU-based acceleration information) at time k. The second term is K (the so-called Kalman gain that minimizes the error covariance of the position at time k + 1) times the difference (residual) between the measured position yk+1 and the prediction of the
426 H.-E. Andersen
measured position (Cxˆ k ). The Kalman gain (K) incorporates the measurement error such that when the measurement error (i.e. GPS error in our case) is large, the gain K is small, and the measured GPS position (yk+1) will not have as much influence on the estimated position xˆ k+1 [50].
Although differential post-processing of the GPS data previously required base station data, even these requirements are now disappearing with the advent of processing techniques such as precise point positioning (PPP), which can provide accurate GPS coordinates without base station data [37]. These developments could dramatically reduce the cost of aerial photo acquisitions, especially in remote, unpopulated areas.
In this section, we first consider manual and then automatic methods of forest photogrammetry.
10.2.2.1Manual Forest Measurements Using Large-Scale Aerial Photogrammetry
The ability to accurately identify and measure individual tree crowns using aerial photographs is heavily dependent on scale and image quality [23]. Scale is a function of focal length of the camera lens and the flying height, while image quality is determined by many factors, including film characteristics and processing, camera lens design, atmospheric and lighting conditions, image motion compensation (blur), resolution (pixel size), color balance, etc.
The other important consideration in determining the efficacy of tree height measurement using forest photogrammetry is image geometry and radial displacement. Many newer digital imaging systems have significantly shorter effective focal lengths than older film mapping cameras, leading to severe layover of trees throughout most of the image area (see Fig. 10.3). Layover (radial displacement) can make it very difficult to view forested scenes in stereo, leading to difficulties in tree height measurement. For this reason, relatively long lenses (e.g. 305 mm) were usually used on mapping cameras in forest photogrammetry applications [23].
Image motion, or blurring due to movement of the aircraft during the time that the camera shutter is open, is another significant concern when acquiring largescale photography for forestry applications. Before the advent of Forward Motion Compensation (FMC) technology in the late 1980s, image motion with film cameras could only be minimized through the use of a short exposure time, which in turn required a fast film with larger grain size [39]. FMC technology actually involved moving the film plate a minute amount during the time that the shutter is open, thereby reducing image blurring. In the case of modern digital imaging systems, a technology called Time Delayed Integration (TDI) is used. In TDI, the image is