Forest canopy height and carbon estimation at Monks Wood National Nature Reserve, UK, using dual-wavelength SAR interferometry

Remote Sensing of Environment 108 (2007) 224 – 239 www.elsevier.com/locate/rse

Forest canopy height and carbon estimation at Monks Wood National Nature Reserve, UK, using dual-wavelength SAR interferometry H. Balzter a,⁎, C.S. Rowland b , P. Saich c a

University of Leicester, Department of Geography, Climate and Land Surface Systems Interaction Center (CLASSIC), University Road, Leicester LE1 7RH, UK b Centre for Ecology and Hydrology (CEH), Section for Earth Observation, Monks Wood, UK c Department of Geography, University College London, 26 Bedford Way, London WC1 0AP, UK Received 27 October 2005; received in revised form 31 October 2006; accepted 4 November 2006

Abstract Forest canopy height is a critical parameter in better quantifying the terrestrial carbon cycle. It can be used to estimate aboveground biomass and carbon pools stored in the vegetation, and predict timber yield for forest management. Polarimetric SAR interferometry (PolInSAR) uses polarimetric separation of scattering phase centers derived from interferometry to estimate canopy height. A limitation of PolInSAR is that it relies on sufficient scattering phase center separation at each pixel to be able to derive accurate forest canopy height estimates. The effect of wavelengthdependent penetration depth into the canopy is known to be strong, and could potentially lead to a better height separation than relying on polarization combinations at one wavelength alone. Here we present a new method for canopy height mapping using dual-wavelength SAR interferometry (InSAR) at X- and L-band. The method is based on the scattering phase center separation at different wavelengths. It involves the generation of a smoothed interpolated terrain elevation model underneath the forest canopy from repeat-pass L-band InSAR data. The terrain model is then used to remove the terrain component from the single-pass X-band interferometric surface height to estimate forest canopy height. The ability of L-band to map terrain height under vegetation relies on sufficient spatial heterogeneity of the density of scattering elements that scatter L-band electromagnetic waves within each resolution cell. The method is demonstrated with airborne X-band VV polarized single-pass and L-band HH polarized repeat-pass SAR interferometry using data acquired by the E-SAR sensor over Monks Wood National Nature Reserve, UK. This is one of the first radar studies of a semi-natural deciduous woodland that exhibits considerable spatial heterogeneity of vegetation type and density. The canopy height model is validated using airborne imaging LIDAR data acquired by the Environment Agency. The rmse of the LIDAR canopy height estimates compared to theodolite data is 2.15 m (relative error 17.6%). The rmse of the dual-wavelength InSAR-derived canopy height model compared to LIDAR is 3.49 m (relative error 28.5%). From the canopy height maps carbon pools are estimated using allometric equations. The results are compared to a field survey of carbon pools and rmse values are presented. The dual-wavelength InSAR method could potentially be delivered from a spaceborne constellation similar to the TerraSAR system. © 2006 Elsevier Inc. All rights reserved. Keywords: Synthetic Aperture Radar (SAR); Interferometry; Tree height; Canopy height; Carbon; LIDAR

1. Introduction The long-term net carbon flux between terrestrial ecosystems and the atmosphere has been dominated by changes in forest area and changes in forest biomass per hectare resulting from management and regrowth (Houghton, 2005). These processes are thus vital for understanding the global carbon balance. The quantification of carbon content in forest ecosystems is important for forest management, climate science and verification of interna⁎ Corresponding author. Tel.: +44 116 252 3820; fax: +44 116 252 3854. E-mail address: [email protected] (H. Balzter). 0034-4257/$ - see front matter © 2006 Elsevier Inc. All rights reserved. doi:10.1016/j.rse.2006.11.014

tional environmental treaties. Forests and woodlands interact with the carbon cycle through photosynthesis, respiration and disturbance events (fire, windfall, insect damage and logging). There are still considerable quantitative uncertainties in the magnitude of the carbon sink in different regions and the contribution of different processes to the overall carbon cycle (Schimel et al., 2001), which could be reduced through better estimates of terrestrial, biotic carbon pools such as forest ecosystems. Increasing awareness of risks associated with climate change has led to international measures trying to limit the increase of greenhouse gas (GHG) concentrations in the atmosphere. The UN Framework Convention on Climate Change and the

H. Balzter et al. / Remote Sensing of Environment 108 (2007) 224–239

Kyoto Protocol aim at a reduction of greenhouse gas emissions to prevent dangerous climate change. Af-, Re- and Deforestation measures are accountable under the Protocol. The possible uses of remote sensing in verifying the Kyoto Protocol are discussed by Rosenqvist et al. (2003). Terrestrial vegetation, mainly forest, is being discussed as a potential carbon sink with the important function of stabilizing GHG concentrations. The Marrakech Accords to the Kyoto Protocol define “forest” as a “minimum area of land of 0.05–1.0 ha with tree crown cover (or equivalent stocking level) of more than 10–30% with trees with the potential to reach a minimum height of 2–5 m at maturity in situ…”. From 2008, nations are required to report human activities that affect their national GHG inventories. To achieve a full carbon account, information is needed on the carbon stocks in forest. In forestry, canopy height is often used to estimate forest biomass and carbon pools, as the quantities are functionally related. Field-based top height measurements are used in estimating timber yield from yield class models by the UK Forestry Commission (Forestry Commission, 1981). The aboveground carbon content of the forest can be estimated using knowledge of forest canopy height and allometric equations (Patenaude et al., 2002). Here, we define forest canopy height as the height of the highest vegetation components above ground level, mean stand height as the mean canopy height of a forest stand, tree height as the height of the tip of the stem of an individual tree, top height of a forest stand as the average total height of the 100 trees of largest diameter at breast height per hectare and the interferometric scattering phase center height from InSAR as the vertical location within the canopy from which most of the backscatter signal is returned. Forest canopy height can be estimated using remote sensing techniques, particularly stereophotogrammetry, LIDAR and Synthetic Aperture Radar (SAR) interferometry (InSAR), including polarimetric interferometry (PolInSAR). Patenaude et al. (2005) give an overview of approaches for above-ground forest carbon stock estimation including RADAR and LIDAR (Light Detection and Ranging). Three basic approaches of biomass mapping from SAR can be distinguished and are described below: approaches based on backscatter, coherence and phase. 1.1. Backscatter based approaches Backscatter intensity is the energy that is received by the SAR sensor after transmission of an energy pulse to the target. The total backscatter from a forest target consists of contributions from a number of basic physical scattering mechanisms, specifically volume scattering from the canopy, rough surface scattering from the ground and double-bounce scattering from trunk-ground interactions. These scattering mechanisms contribute to the return signal at specific polarization combinations, and depend amongst other factors on the electromagnetic wavelength. Radar backscatter intensity has been used to estimate woody biomass of forest (Kasischke et al., 1997), forest biomass (Le Toan et al., 1992), forest aboveground dry biomass (Rignot et al., 1994) and timber volume (Balzter et al., 2002a). Backscatter typically increases with increasing forest biomass, but

225

this function saturates at a wavelength dependent biomass density, which limits the usefulness of SAR for this application (Imhoff, 1995). The form of the functional relationship between backscatter and biomass depends heavily on vegetation structure, which can confuse the retrieval of biomass. A backscatter modeling study using MIMICS-I for L-band radar for dry and wet peat swamp forest found that the canopy layer contributed significantly to radar backscatter, but some L-band radiation also penetrated through the canopy to the ground layer, as shown by significant trunk-ground contributions in both flooded and nonflooded peat-swamp forest (Aziz & White, 2003). A multitemporal analysis of SEASAT and JERS-1 SAR over a coniferous forest plantation showed that L-band backscatter difference is sensitive to tree growth over an 18 year period (Balzter et al., 2003). Hyyppä and Hallikainen (1996) present results from the helicopter-borne HUTSCAT ranging scatterometer at C- and X-band (polarimetric) over Finnish boreal forest stands, which showed a tree species dependent retrieval accuracy of tree height. Tree height of pine stands could be estimated more accurately than deciduous, spruce or mixed stands. From 3° to 40° incidence angle the retrieval accuracy decreased. The profiling capabilities of this instrument mean that ground elevation can be estimated from the signal, even though a certain amount of noise means that digitizing the surface can be more accurate. Martinez et al. (2000) relate HUTSCAT data to a radiative transfer model coupled to a tree structure model. 1.2. Coherence based approaches Coherence based approaches are based on the estimation of the complex correlation coefficient between two SAR acquisitions. The coherence is the magnitude of the complex correlation coefficient, while the interferometric phase is defined as an angular measurement. Coherence based approaches are based on the observation that interferometric coherence under certain conditions is related to forest biomass. Forest stem volume was retrieved from multitemporal interferometric ERS 1 and 2 C-band tandem coherence under winter conditions over forest test sites in Finland and Sweden by Askne and Santoro (2005), who conclude that the error level for large forest stands (N 2 ha) in a well managed and homogeneous boreal forest may be expected to be in the 15% to 25% range, deteriorating for small and heterogeneous stands. The superior performance of winter coherence over summer coherence for stem volume retrieval was also found by Pulliainen et al. (2003). At L-band, multitemporal JERS-1 data showed that interferometric wintertime coherence with a 44 day temporal baseline is suitable for estimations of growing-stock volume in boreal forest in Siberia (Eriksson et al., 2003). Tansey et al. (2004) present a forest growing stock classification algorithm based on the information content of both ERS tandem coherence and JERS backscatter that gives classification accuracies of greater than 70% for test sites in Sweden, Siberia and the UK. 1.3. Phase based approaches Phase-based InSAR techniques exploit the interference patterns (so-called fringes) of two electromagnetic waves to

226

H. Balzter et al. / Remote Sensing of Environment 108 (2007) 224–239

remove the random component of the phase of the radar signal to estimate the topographic height of the scattering phase center. The location of the scattering phase center within the canopy depends on vegetation structure, scattering mechanisms and sensor characteristics. Over forested areas, the scattering phase center is the integral of the returns from a large collection of scatterers, which include ground, stems, branches and leaves or needles. Which types of scatterers interact most strongly with the radar wave depends on the wavelength, polarization, incidence angle and vegetation density. Interferometric techniques can be used for estimating forest canopy height because the interferometric phase relates to the scattering phase center, which at suitable wavelengths and polarizations is composed of terrain height and vegetation canopy height (Balzter, 2001). Phase unwrapping is an important step in this processing chain, and a wide range of methods have been developed to get good phase estimates (Bamler et al., 1998; Costantini, 1998; Costantini et al., 1999; Fornaro & Sansosti, 1999; Fornaro et al., 1996, 1997; Griffiths & Wilkinson, 1994; Lin et al., 1994; Lyuboshenko & Maitre, 1999; Pritt, 1996; Trouve et al., 1998; Wang & Li, 1999; Xu & Cumming, 1999; Zebker & Lu, 1998). If there are significant gaps in the canopy, the interferometric height of the scattering phase center will be a composite of the canopy top and the proportion of the ground that contributes to the signal in the resolution cell (Hagberg et al., 1995). The use of InSAR for estimating canopy height (Balzter, 2001) was demonstrated using C-band interferometric height discontinuities at forest edges by Hagberg et al. (1995) and using effective interferometric C-band height from ERS-1 by Askne et al. (1997). Treuhaft et al. (1996) carried out a modeling study for airborne C-band data and found that for an accuracy of 0.5% in the InSAR normalized cross-correlation amplitude and an accuracy of about 5° in the interferometric phase, vegetation layer depths of a few meters and ground surface heights can be estimated from InSAR for many types of vegetation layers. Treuhaft et al. (1996) found an average agreement of InSAR derived vegetation layer depth with ground data of 5 m, with the error being related to stand height. Interferometric SIR-C shuttle radar data at C- and L-band were studied over tropical rain forest by Rignot (1996). At L-band, the rms difference in inferred topographic height between the forest and adjacent clearings was 5 m, equivalent to the height noise (Rignot, 1996). The dual-wavelength SAR system TOPOSAR is an airborne sensor system consisting of a single-pass X-band and a repeat-pass P-band interferometric SAR constellation. Andersen et al. (2004) reported that this system can be used to estimate mean stand canopy height, canopy fuel weight and other biophysical parameters. Repeat-pass InSAR is limited in its accuracy by temporal decorrelation caused by temperature, orientation angle and moisture changes between the two SAR acquisitions (Hagberg et al., 1995; Zebker & Villasenor, 1992), and depends on the wavelength (Balzter, 2001; Zebker & Villasenor, 1992). A study conducted by Kellndorfer et al. (2004) to determine the feasibility of vegetation canopy height retrieval from digital elevation data collected during the 2000 Shuttle Radar Topography Mission (SRTM) with a C-band InSAR instrument suggested that a minimum mapping unit of

approximately 1.8 ha is needed to achieve good estimates. Approaches have been developed that take into account forest structure in biomass estimation. Forest biomass was estimated from multi-altitude, airborne, C-band SAR interferometry density profiles, which were normalized by leaf area index from airborne hyperspectral optical imagery, yielding the forest canopy leaf area density (Treuhaft et al., 2003, 2004). One approach to estimate forest canopy height is using fully polarimetric interferometry (PolInSAR). PolInSAR methods (Cloude & Papathanassiou, 1998) make use of full polarimetry to estimate components of different scattering mechanisms. After applying an eigenvector-based coherence optimization algorithm (Cloude & Papathanassiou, 1998) the scattering phase centers of the three optimized channels show better separation in the complex unitary circle compared to the HV (horizontal transmit, vertical receive) basis (Stebler et al., 2002). A coherent scattering model for single-baseline polarimetric SAR interferometry (PolInSAR) is described in (Papathanassiou & Cloude, 2001). The inverted model can be used to estimate canopy height, average extinction and terrain elevation. A limitation of PolInSAR is that it relies on sufficient scattering phase center separation at each pixel to be able to derive accurate canopy height estimates. The effect of wavelengthdependent penetration depth into the canopy is known to be strong, and could potentially lead to a better height separation than relying on polarization combinations at one wavelength alone. In this paper we describe the derivation of a spatial map of forest carbon pools based on a canopy height model derived through a number of processing steps from dual-wavelength SAR interferometry. This technique is motivated by coherent microwave modeling results, which are presented in Section 3.1, indicating that for L-band HH polarization the interferometric scattering phase center height can provide an estimate of terrain elevation under certain conditions (depending on the density and size of scattering elements in a resolution cell, ground properties, viewing geometry and moisture conditions) and the signal at Xband VV polarization if acquired in single-pass interferometric mode can be used with certain adjustments made for penetration into the canopy as a basis to estimate mean canopy height. 2. Data and methods 2.1. Remote sensing data and processing methods 2.1.1. Test site description Monks Wood NNR is a long-term calibration and validation site for the Airborne Research and Survey Facility operated by the UK Natural Environment Research Council. It is also used for long-term monitoring of plant, butterfly, moth and bird populations. The reserve covers 157 ha of mostly temperate deciduous woodland and is located at 52°24′ N, 0°14′ E in Cambridgeshire. The highest part of the reserve rises to 46 m above mean sea level with a maximum slope angle of 14.5°. The dominant tree species are ash (Fraxinus excelsior), oak (Quercus robur) and field maple (Acer campestre) with some occasional clusters of elm (Ulmus carpinifolia) and aspen

H. Balzter et al. / Remote Sensing of Environment 108 (2007) 224–239

(Populus tremula) on wetter soils. Dominant shrub species include hawthorn (Crataegus monogyna), hazel (Corylus avellana), blackthorn (Prunus spinosa), dogwood (Cornus sanguinea), and wild privet (Ligustrum vulgare). The vegetation canopy is diverse and species composition, canopy density and height vary substantially across the extent of the nature reserve resulting in high spatial heterogeneity in scattering properties. The wood is not subject to thinning practices but is managed for nature conservation. Strong westerly winds prevailed on the acquisition date, and the flight path had to be adjusted to fly in the wind direction to be able to maintain an accurate baseline for the repeat-pass interferometry. The total vegetation biomass carbon content varies between about 80 and 210 t/ha with an additional 270 to 400 t/ha carbon stored in the soil (see Patenaude et al., 2003 for details). 2.1.2. SAR dataset and processing Airborne SAR data were acquired over Monks Wood National Nature Reserve (NNR), East Anglia, UK, during the SAR and Hyperspectral Airborne Campaign (SHAC) on 1 June 2000 by the E-SAR sensor operated by the German Aerospace Research Institute (DLR). The E-SAR was flown in repeat-pass mode along parallel tracks at a baseline of about 10 m. L-band was acquired in quad-polarized mode (repeat-pass interferometry) and X-band in VV polarized mode (single-pass interferometry with 1.5 m baseline). Table 1 gives details of the acquisition parameters. The complex correlation coefficient at each pixel was derived by forming the complex conjugate of the two SAR acquisitions. The absolute interferometric phase induced by the terrain for the SAR viewing geometry is (Madsen and Zebker, 1998): /¼

2p pDR k

ð1Þ

where p = 1 when the channels share the transmit antenna, and p = 2 when each channel is transmitting and receiving using its

Table 1 E-SAR acquisition parameters for L-band (top) and X-band (bottom) Longitude Easting (Start of flight track) Latitude Northing (Start of flight track) Acquisition date

− 0.28755888° 52.429474° 01-JUNE-2000

Band

L

Frequency Polarisation InSAR Mode Baseline length Range pixel spacing (SLC) Azimuth pixel spacing (SLC)

1.3 GHz quad-pol repeat-pass ∼ 10.0 m 1.49854 m 0.853503 m

Band

X

Frequency Polarisation InSAR Mode Baseline length Range Pixel Spacing (SLC) Azimuth Pixel Spacing (SLC)

9.6 GHz VV single-pass 1.5 m 1.49854 m 0.357801 m

227

own antenna. Due to decorrelation effects the measured interferometric phase ϕM is disturbed by noise. The phase noise ϕN may be considered as being additive to the measured phase ϕM and thus cause it be wrapped in the interval − π to π (Bamler & Hartl, 1998). For terrain reconstruction from acrosstrack interferograms this ambiguity must be resolved during the phase unwrapping (Madsen and Zebker, 1998). Let the operator W{.} that produces the wrapped phase ϕM be defined as: /M ¼ W f/g ¼ modf/ þ p; 2pg−pað−p; pÞ

ð2Þ

The aim of phase unwrapping is to find an estimate ϕˆ of the ‘true’ phase ϕ given its principal wrapped value ϕM. The phase gradient ∇ϕM is the quantity where disturbing contributions from the wrapping operator of equation 2 can be separated from the true phase. This argument leads to the following two-step paradigm underlying the majority of phase unwrapping algoˆϕM of the phase gradient is obtained rithms: (i) an estimate ∇ from the wrapped phase; (ii) this estimate is integrated either via a one-dimensional summation or by a 2D convolution where weighting according to the reliability of gradient estimates is applied in some cases. Phase unwrapping algorithms have different approaches how to compute these two steps (Bamler & Hartl, 1998). The interferometric correlation is a measure of the accuracy of the estimation of the interferometric phase. Interference phenomena such as fringes will be observed as long as there is at least partial coherence between the two images. With very low coherence, the phase is not derived accurately, and an error in phase translates directly into an error in height measurement. γ can also be expressed as the product of the three components of system noise and the differences in target properties over the temporal and spatial baseline (Zebker & Villasenor, 1992): jgj ¼ jgnoise jd jgbaseline jd jgtemporal j

ð3Þ

|γnoise| is related to the signal-to-noise ratio (SNR) of the sensor, and is only significantly different from unity for areas of very low backscatter (e.g. open water bodies): jgnoise j ¼

1 1 1 þ SNR

ð4Þ

|γbaseline|, the coherence term representing baseline decorrelation, arises from surface and volume scattering and can be separated as follows: jgbaseline j ¼ jgslant

range jd

jgvolume j

ð5Þ

|γslant range| can be increased to 1 by applying common band filtering, as was done in the present study. However, |γvolume| is close to 1 only in the case when the scatterers are confined to a plane. When they are distributed in depth, as for example with multiple scattering from a forest canopy, then there is appreciable volume decorrelation (Balzter, 2001). |γtemporal| is the coherence term representing decorrelation of the target arising from the time separation of the two observations e.g. when the images have been taken using repeat-pass system. It is close to unity for stable, human-made structures, moderately reduced for

228

H. Balzter et al. / Remote Sensing of Environment 108 (2007) 224–239

bare soil surfaces and agricultural crops and can be substantial for tall vegetation such as forests. In the latter case, temporal changes are caused by changes of the scatterers (growth or loss of foliage, wind motion) and changes in their dielectric constant (surface moisture, freezing, thawing). From Eq. (3) it can be deducted that because |γnoise| is only significantly different from unity for areas of very low backscatter and could be assumed 1 for the test site described in this study: jgbaseline j ¼ jgslant

range jd

jgvolume j

ð6Þ

|γslant range| has been increased to 1 by applying common band filtering, thus it follows as an approximation that jgj ¼ jgvolume jdjgtemporal j

ð7Þ

Two approaches can be used to estimate the volume and temporal decorrelation components of the observed coherence, both approaches are model-based. (i) The Interferometric Water Cloud Model can be used to model volume decorrelation. For large two-way attenuation and when the normal component of the baseline is below 100 m, the volume decorrelation at C-band becomes negligible in this model (Santoro et al., 2003), so in actual fact the observed coherence then reflects primarily temporal decorrelation. (ii) The volume decorrelation can be estimated using the method of Gaveau (2002) who proposed a simple modelbased methodology for estimating the level of volume decorrelation at C-band, for a given canopy, and hence indirectly estimating the temporal decorrelation. Both approaches were developed for C-band, and make assumptions about the vertical distribution of scattering in the canopy which is only sufficiently quantified at C- and X-band but not yet at L-band. First we have to define three different types of elevation models: Digital Surface Models (DSM) are defined as the elevation surface from the remote sensing data which include a terrain and a vegetation term, Digital Terrain Models (DTM) are defined as the elevation of the unvegetated terrain, and Canopy Height Models (CHM) as the elevation relating to radiation processes in the vegetation canopy. DSM's are closely related to the scattering processes that formed the SAR images and the reflection processes that underlie the LIDAR image, DTM's may be derived by interpolation of a selection of pixels based on the vertical proximity of the scattering to the underlying terrain and can contain interpolated and filtered data, and CHM's are generally based on the subtraction of a DTM from a DSM with some adjustment for penetration depth that is vegetation type dependent. In this study, firstly DSMs were derived from the X-band VV polarized single-pass interferometric SAR data, and from the Lband HH polarized repeat-pass interferometric SAR data with a 13 min temporal gap between acquisitions (11:21am and 11:34 am). The interferometric processing involved common band

filtering, interferogram generation with multi-looking (2 looks in range and 6 looks in azimuth direction), curved Earth phase trend removal, and coherence estimation with an adaptive window size between 2 and 5 and two possible weighting functions. The adaptive coherence estimation algorithm in the Gamma Interferometric SAR Processor (Wegmüller & Werner, 1997) improves the estimation accuracy of low coherence values by increasing the estimation window for pixels with low coherence and thus reducing the bias. For the coherence estimation of each pixel a corresponding window size and weighting function was selected. Then within the determined window the coherence of the central pixel was estimated as a weighted average of all pixels with local weights being related to the magnitude of the pixelwise coherence. In the case of an initial high coherence estimate of the central pixel, the coherence was estimated with a small window and areas of a high initial estimate should only contribute to the estimation process at a short distance. In the case of an initial low coherence estimate, a large window was applied. Following the coherence estimation, a phase unwrapping algorithm was applied. The flattened phase was filtered adaptively to the local slope using a Kaiser window function with a width of 6 pixels and a signal-tonoise ratio threshold of 0.2 to reduce the phase noise. Areas with coherence less than 0.25 were masked out. Phase residues were determined to avoid inconsistencies in the unwrapped phase image. Phase unwrapping was then performed using a region growing algorithm starting in the image center. To derive a digital elevation model unwrapped phase values were extracted for manually defined vegetation-free ground control points with known topographic height. The interferometric baseline was then estimated for each image line using a least-squares procedure, which also provided the rmse of the DSM at the ground control points. The unwrapped phase is a measure of the propagation path length difference between the two radar pulses. It was combined with the baseline data for each image line to calculate the look vector from the radar to the point on the ground. Terrain height was estimated based on this acquisition geometry information. A transformation from slant range (SAR acquisition geometry) to ground range (pixels have equal area on the ground) was the final step of the processing chain. An iterative resampling was carried out with window sizes of 7, 15 and 31 pixels to fill gaps in the ground range interferometric scattering phase center height images. All DSMs were coregistered and resampled to 3 m pixel spacing. 2.1.3. LIDAR An airborne imaging LIDAR acquisition by the Environment Agency is used for assessing the accuracy of the SAR-derived products. The Airborne Laser Terrain Mapper (Optech ALTM 1210) was flown on the 10th of June 2000. On average, one laser point was recorded every 4.83 m2. The LIDAR spot heights have an accuracy of about 15 cm according to technical specifications, which was verified in a comparison with Total Station measurements of unvegetated areas at Monks Wood (see Hill et al., 2002, for details). The retrieval accuracy of an interpolated DTM under forest from LIDAR last return heights at the Monks Wood test site has been published by Hill and

H. Balzter et al. / Remote Sensing of Environment 108 (2007) 224–239

Thomson (2005) and is approximately 0.5 m. The first and last returns of the laser waveform were recorded for an irregularly spaced grid of points. From the laser point height measurements, a gridded digital CHM and a DTM of the underlying ground were generated from the first and last return respectively as described in Patenaude et al. (2004). The LIDAR images were coregistered with the SAR images and resampled to 3 m pixel spacing. The mean, median, mode and maximum canopy height of the wooded area at Monks Wood derived from the LIDAR CHM excluding pixels below 2 m are 12.23 m, 12.98 m, 14.5 m and 26.29 m respectively. 2.2. Estimating forest canopy height and carbon pools The basic principle of directly estimating canopy height is to subtract the terrain elevation from the surface elevation retrieved from remote sensing. The availability of a good terrain model is often a limitation of this approach. To generate the dual-wavelength InSAR CHM, (i) the X-band InSAR scattering phase center heights provide an initial digital surface model (DSM), (ii) a smoothed and interpolated digital terrain model (DTM) is generated from a subset of pixels with sufficiently low density of scattering elements in the resolution cell to provide a signal from near the terrain surface at the L-band HH polarization, and (iii) the terrain model is subtracted from the surface model to provide a canopy height model (CHM). The algorithm to generate the DTM underneath the forest canopy from the L-HH interferogram is as follows: (i) Select all pixels with L-HH coherence greater than 0.6 in the L-HH InSAR DSM. These pixels are assumed to have a strong ground or trunk/ground interaction contribution to the radar signal as a result of sufficiently low density of vegetation components scattering the radiation at L-band. The temporal decorrelation caused by the movement of scatterers between the two acquisitions caused a loss of coherence over densely vegetated areas as a result of a stronger displacement of these vegetation components between the repeat-pass acquisitions and volume decorrelation, which allows identification of low scattering element density pixels. This step results in a selection of highcoherence pixels. (ii) Interpolate the interferometric height estimates in the DSM for the selected high-coherence pixels using a radial search algorithm to fill in masked-out lowcoherence pixels. For each low-coherence pixel, a search radius is used to find high-coherence DSM values within a circular neighborhood. If the number of pixels that can be used for the interpolation is greater than or equal to a specified minimum number, the unweighted average of all points within the circle is calculated and assigned to the pixel. If not, the radius is extended up to a set limit. Here the algorithm started with a search radius of 10 pixels and stepwise increased the radius to a maximum of 40 pixels until at least 30 high-coherence DSM values were found that could be used to provide the interpolated value of the low-coherence center pixel. (iii) Smooth the interpolated elevation data product using a 21 by 21 moving average filter to reduce noise. The smoothed interpolated terrain elevation model will be called the DTM.

229

This terrain retrieval algorithm is based on the physical scattering processes in forest canopies which exhibit a sufficient degree of spatial heterogeneity of the density of scattering elements. Electromagnetic waves are predominantly scattered in forest canopies through a wavelength and polarization dependent mix of the three scattering mechanisms of volume scattering, double-bounce scattering and surface scattering. The increased temporal decorrelation over forest due to the movements of scatterers (branches, leaves and twigs) between the two L-band acquisitions is suggested as one of the reasons why the selection of high-coherence pixels mainly includes spatial areas in which the scattering phase centers were relatively close to the terrain underneath the (sparse) forest canopy. Coherence of the double-bounce and surface scattering mechanisms is likely to be higher than for volume scattering, because stems and thick branches do not move as much as twigs and leaves in a short time interval, and the ground roughness and moisture may be assumed stable over a 15 min time period. This suggests that in the case of the Monks Wood E-SAR data high-coherence pixels will tend to represent high proportions of double-bounce and ground returns and low proportions of canopy returns. The coherence threshold was experimentally determined to be the highest coherence threshold that still yields enough pixel values to avoid an increasing error from interpolating over too large masked out low-coherence areas. The validated DTM from the LIDAR image was used for comparison. Its accuracy was determined to be 0.5 m (Hill et al., 2002). When the rmse between a set of variants of the smoothed interpolated L-HH DTM and the LIDAR DTM is calculated for a number of coherence thresholds between 0.0 and 0.9 it becomes apparent that the rmse does not change much for coherence thresholds between 0.0 and 0.6, but increases rapidly for coherence larger than 0.6. This behavior may be explained by the fact that on one hand, for low coherence thresholds, many pixels are included in the interpolation but they have higher errors when compared to the ‘true’ terrain values. On the other hand, for high coherence thresholds a small number of pixels are retained for the interpolation but the terrain estimates are more accurate at those pixels. So by increasing the coherence threshold the error component from interpolating over masked out pixels increases, whilst at the same time the error component introduced by high phase noise and stronger canopy contributions to the signal at low coherence pixels decreases. By the time the coherence threshold increases above 0.6 there are insufficient highcoherence pixels retained for interpolation in this dataset to produce a surface which resembles the terrain. Below a 0.6 threshold the retained LHH DTM pixels are spatially correlated, so there is a degree of spatial redundancy, and reducing the coherence threshold further has minimal impact on the derived DTM. When comparing different variants of the LHH DTM for coherence thresholds below 0.6 it also becomes apparent that although they show similar overall rmse values, the spatial distribution of errors changes for different thresholds. The CHM from dual-wavelength SAR interferometry was estimated as the difference between the X-band DSM and the smoothed interpolated L-band DTM. The dual-wavelength InSAR CHM was adjusted with offset constants of 1.4m

230

H. Balzter et al. / Remote Sensing of Environment 108 (2007) 224–239

LIDAR CHM based on a subtraction of the LIDAR DTM from the X-band DSM and a second a LIDAR CHM based on a subtraction of the LIDAR DTM from the LIDAR DSM. The carbon estimation follows allometric methods applied to the remotely sensed canopy height models. All constants and factors mentioned in this section are described in detail in Patenaude et al. (2002). The CHM's from InSAR and LIDAR were used to derive top height for all trees taller than 6 m by applying an allometric adjustment factor of 1.0492 based on published tables by Hamilton (1975). Stem volume was derived from top height using Eq. (3) from Patenaude et al. (2002) which is based on general yield class models used by the UK Forestry Commission: Vstem ¼

6:357834 expð0:174045dhtop Þ 0:59822

ð8Þ

where Vstem is the stem volume (m3/ha) and htop is the top height of the stand (m). This equation assumes a constant basal area/top height ratio.

Fig. 1. Scattering phase center height at X-band VV polarization and L-band HH, HV and VV polarization (H=horizontal, V=vertical). (a) from the coherent microwave model CASM for tree densities of 1000 (dotted lines) and 2000 ha−1 (solid lines) at 45° incidence angle. The maximum phase center separation is between the phase centers of X-VV and L-HH. Reproduced in a slightly modified form from Balzter et al. (2002b). (b) from E-SAR data over Monks Wood National Nature Reserve (mean ± standard deviation is shown for X-VV and L-HH). The mean scattering phase center heights in (b) were estimated from the E-SAR InSAR X- and L-band data by subtracting the LIDAR DTM from the InSAR DSM's and averaging over all pixels within the respective height class. The data were resampled from 3 m × 3 m to 12 m × 12 m. The number of pixels averaged varies from 13 for the 22 m class to1023 for the 15 m class.

(CHM ≤ 4 m) or 2.25 m (CHM N 4 m) that were derived for the LIDAR derived canopy height model by Patenaude et al. (2002) to adjust for penetration into the canopy through gaps in the foliage. As Fig. 2 demonstrates, the X-VV scattering phase center heights in the interferometric DSM are of similar magnitude to the LIDAR first return heights. The infrared pulse from the LIDAR that is recorded as the first return signal is reflected from the first object it encounters, which is usually near the upper edge of the canopy if the canopy is sufficiently closed. Similarly the X-band SAR signal is scattered by objects approximately equivalent to the wavelength, such as leaves and small twigs. The vertical heights of the reflection and scattering processes are thus very similar. Because of this observation, the abovementioned adjustment offsets determined from LIDAR were assumed to be valid for the X-band InSAR CHM as well. For comparison two more CHM's were derived: an X-band/

Fig. 2. Profile plots of elevation data along the transect indicated in the top left corner. (a) digital surface models (DSM), (b) canopy height models (CHM) derived by subtracting the LIDAR last return digital terrain model (DTM) from the DSM's in (a). The CHM's here are models of the vertical location of the scattering processes that influence the remote sensing signals. The scattering phase centers at XVV are close to the LIDAR reflecting elements at the canopy top, while L-band penetrates deeper into the canopy at all polarizations.

H. Balzter et al. / Remote Sensing of Environment 108 (2007) 224–239

Carbon content was estimated by multiplying stem volume with a biomass expansion factor for deciduous British forest of 1.36, a dry mass conversion factor of 0.55 and a carbon content conversion factor of 0.49 from Patenaude et al. (2002). The ratio between overstorey and understorey carbon content is a negative exponential function of overstorey carbon content and can thus be estimated (Patenaude et al., 2002). For canopy height greater than 15 m the understorey carbon content was estimated as follows:

231

The sum of both carbon pools is defined as the total vegetation carbon pool within a spatial resolution cell (pixel). The carbon maps were averaged over a 7 × 7 window to 21 m pixel spacing, which meets the requirements of 0.05 ha minimum reporting unit outlined in the Marrakech Accord to the Kyoto Protocol. 3. Results 3.1. Estimating terrain elevation from L-band HH InSAR

Cunderst ¼ Coverst d ð59:133d C−1:2977 overst Þ

ð9Þ

where Cunderst and Coverst are the carbon content (t/ha) of the understorey and overstorey respectively.

The coherent microwave scattering model CASM (Saich et al., 2001) is based on the coherent addition of single scattering from vegetation elements, along with double-bounce

Fig. 3. Height maps of Monks Wood National Nature Reserve. (a) L-HH InSAR DSM from scattering phase center heights; (b) L-HH InSAR DSM pixels with coherence b0.6 masked out (black); (c) interpolated and 21 × 21 average-filtered DTM based on (b); (d) interpolated LIDAR DTM. Color scale: 0…70 m. Area about 1 × 1 km.

232

H. Balzter et al. / Remote Sensing of Environment 108 (2007) 224–239

Fig. 4. Validation of the smoothed interpolated L-HH Digital Terrain Model (DTM) by comparison with theodolite measurements for open areas (‘o’ symbols, e.g. bridleways without tree cover), and for locations underneath the canopy (‘+’ symbols).

scattering between the vegetation and ground and a ground term. The scattering amplitudes are weighted by a term incorporating the phase change due to the propagation path for all the elements in a simulation cell and are then summed to calculate the backscatter intensity. Backscattering and interferometric signatures are obtained by averaging over a large number of independently generated simulation cells. As such, the model allows for prediction of both the backscattering coefficient and complex interferometric coherence. A modeling experiment with CASM indicated that among the polarizations and wavelengths that were acquired during the SHAC campaign, L-band

HH polarization can, in certain conditions, penetrate sufficiently deep into the canopy to provide estimates of terrain elevation under sparse forest canopies, with sufficiently low density of scattering elements active at L-band. Radiation at X-band VV polarization is under certain conditions scattered in a layer close to the canopy height (Fig. 1a). CASM represents tree structures as an assembly of individual scatterers (stem, first and higher order branches, leaves or needles). A tree structural model based on a Lindenmayer system (Prusinkiewicz & Lindenmayer, 1990) was used to run CASM, so the model results shown here are only illustrative and not tailored to the deciduous vegetation at Monks Wood. The modeled wavelength-specific phase scattering centers follow a similar pattern as the observed values at Monks Wood NNR (Fig. 1b). Differences in penetration depth between X- and L-band were also found in a transect analysis of the airborne SAR data (Fig. 2). However, the polarization effects on the scattering phase center heights at L-band that were expected from the CASM model were not detectable in the data (Fig. 1). The modeled HV and VV polarizations at L-band in Fig. 1a show higher scattering phase centers than those observed in Fig. 1b. This could be caused by an overrepresentation of vertically aligned scatterers in the Lindenmayer system tree model causing a too strong vertical scattering component in the model, and needs further improvement to become more realistic. In Fig. 1b the mean L-HH signal originates from slightly above the terrain under the forest. Double-bounce contributions to the SAR signal from trunk-ground interactions may be an explanation of this offset. The L-HH returns also exhibit high variability between different pixels, which leads to the suggestion that scattering processes at L-band over the test site are spatially highly heterogeneous, and in conclusion the distribution of scattering element density in the resolution cells

Fig. 5. Digital Surface Models (DSM's) of Monks Wood National Nature Reserve. (a) X-VV InSAR DSM; (b) LIDAR DSM; color scale: 0…70 m above mean sea level.

H. Balzter et al. / Remote Sensing of Environment 108 (2007) 224–239

233

Fig. 6. Canopy height models (CHM's) of Monks Wood National Nature Reserve. (a) dual-wavelength InSAR CHM; (b) LIDAR CHM (validated using theodolite: rmse=2.15 m). Color scale: 0…25 m.

is highly variable. The standard deviations for L-HH in Fig. 1b show that a proportion of pixels within each tree height class originate from scattering processes sufficiently close to the forest floor to attempt an interpolation of an underlying terrain model. This subset of pixels is probably characterized by a sufficiently small density of scattering elements at L-band in the resolution cell. The spatial heterogeneity of vegetation density is assumed to be a prerequisite for applying the described DTM derivation algorithm. Even over the short 13 min temporal baseline, decorrelation effects in the L-band DSM's were visible at all polarizations over forested areas. The herbaceous vegetation (transect positions from 75 m to 186 m in Fig. 2) produced less noisy scattering phase center height estimates at L-band than the forest area (186 m to 408 m on transect), as was to be expected because of the volume decorrelation within the forest canopy. The X-VV CHM in Fig. 2b is very similar to the LIDAR CHM (rmse = 3.49 m). Fig. 3a shows the L-HH interferometric scattering phase center height map for the area of interest in this study. The full SAR image covers a larger area, but was cut to focus on the Nature Reserve area, and to exclude problematic areas with high incidence angles from the phase unwrapping. Fig. 3b shows all pixels with a coherence greater than 0.6 (other pixels are masked out in black), leaving only pixels with a highly coherent signal probably caused by a strong ground or trunk-ground interaction contribution to the overall signal. These pixels were then interpolated and smoothed to produce the L-HH derived DTM presented in Fig. 3c. For comparison, the LIDAR DTM is shown in Fig. 3d. The interpolated L-band DTM is able to identify major topographic features that are also apparent in the more detailed LIDAR DTM. A comparison of the X-VV DSM with the LIDAR DSM resulted in an rmse of 3.26 m. The L-HH

DSM and the LIDAR DSM show a much higher rmse of 11.43 m, which can be explained by the greater penetration into the canopy at L-band compared to X-band. A comparison of the L-HH DTM with the LIDAR DTM gave an rmse of 3.81 m. The validity of the L-HH DTM was also assessed against a set of theodolite derived terrain height measurements (Fig. 4). It is evident that the DTM represents an approximation to the underlying terrain underneath the forest canopy. Some of the slight overestimation of the theodolite-derived terrain height by the L-HH DTM is thought to originate from double-bounce scattering between trunks and the ground, which can affect both the “under canopy” points in Fig. 4 and those “ground” points that are adjacent to trees and where the viewing angle leads to such a double-bounce scattering.

Fig. 7. Variation of the observed rmse of the dual-wavelength InSAR CHM compared to the LIDAR CHM in range direction for Monks Wood NNR. rmse is calculated linewise over all pixels within blocks of 10 adjacent lines with equal range distance, excluding masked areas. The rmse increases with ground range distance (see regression line). The total rmse is indicated as a dashed line (rmsetot = 3.49 m).

234

H. Balzter et al. / Remote Sensing of Environment 108 (2007) 224–239

3.2. Estimating forest canopy height and carbon pools The X-VV DSM is shown in Fig. 5a, with the LIDAR first return heights in Fig. 5b for comparison. The X-band scattering phase center height estimates in the DSM are not affected by temporal decorrelation because they were acquired in a single pass. Fig. 6a shows the dual-wavelength InSAR CHM compared to the LIDAR CHM in Fig. 6b. The dual-wavelength InSAR CHM tends to provide lower estimates than the LIDAR CHM, as is evident from areas of green to yellow pixels in Fig. 6a with corresponding red pixels in Fig. 6b. The observed rmse of the dual-wavelength InSAR CHM compared to the LIDAR CHM as a function of range position is shown in Fig. 7. The rmse for the first 100 m of around 1.5 m in Fig. 7 is lower than the average because this is an agricultural field. After 100 m the rmse shows an increasing trend with increasing range position due to higher phase noise in the far range. The terrain

elevation of unvegetated areas can be estimated more accurately than in the presence of a canopy. The overall rmse of the dualwavelength InSAR CHM for the area shown in Fig. 6 is 3.49 m. This rmse can have a high impact on the carbon accounting because of the error propagation when converting canopy height to stem volume and then carbon content. Using the carbon accounting method described in Patenaude et al. (2002) this height rmse would imply a resulting carbon rmse of 175 tC/ ha for tall stands of around 25 m canopy height but only 8.5 tC/ ha for young stands of 1–3 m height. To test whether the deviation of the dual-wavelength InSAR CHM varies with canopy height, the root mean square difference (rmsd) between this CHM and the LIDAR CHM was computed for a sample of 1683 pixels. The rmsd increases with canopy height, showing a weak (r2 = 0.066) but significant ( p b 0.0001) relationship. A linear regression resulted in rmsd = 1.5 + 0.13 ⁎ LIDARCHM; excluding canopy heights below 3 m.

Fig. 8. Images of residual height between (a) Lidar-CHM- InSAR-CHM, (b) X-DSM minus Lidar DSM, (c) LHH-DTM- Lidar-DTM. The boundaries of the forested area are indicated as a white line.

H. Balzter et al. / Remote Sensing of Environment 108 (2007) 224–239

Maps of the residual height between the LiDAR derived CHM, DSM and DTM compared to the radar derived equivalents are given in Fig. 8. Bridleways and forest edges sometimes appear white (high difference between the LiDAR and radar derived height), indicating effects of layover and shadow in the radar image. Fig. 9 shows maps of carbon content derived from three different approaches to terrain removal from a DSM. The CHM's from the X-VV interferometry DSM and the LIDAR DSM, which both use the LIDAR DTM to subtract topography, produce similar patterns of carbon distribution (Fig. 9a, b). Fig. 9c shows the carbon stocks from the dualwavelength InSAR CHM, which reproduces similar patterns, although particularly in the far range (bottom of Fig. 9c) the estimates get noisier due to increasing incidence angles. To assess the correspondence of the remote sensing derived carbon

235

pool estimates from the LIDAR CHM and the dual-wavelength InSAR CHM, carbon pools were calculated for 5 selected stands at Monks Wood. These five values were compared to carbon pool estimates from a previous field survey of the same stands published by Patenaude et al. (2003) using the rmse. The rmse of the aboveground carbon content per hectare derived from the LIDAR DSM / LIDAR DTM was estimated to be 34.9 tC/ha, while it was 54.5 tC/ha from the X-band DSM / LIDAR DTM, and 69.0 tC/ha from the dual-wavelength InSAR CHM when all 5 stands were included in the calculation. However, the airborne L-band interferogram is strongly affected by phase noise in the far range of the image due to large incidence angles, and stands 1 and 4 fall outside the region for which the interferogram shows good coherence. After eliminating these stands as outliers, the rmse for the dual-wavelength derived carbon content

Fig. 9. Vegetation carbon content [t/ha] at Monks Wood NNR derived from the canopy height models from (a) LIDAR DSM and LIDAR DTM, (b) XVV InSAR DSM and LIDAR DTM, (c) XVV InSAR DSM and smoothed interpolated LHH InSAR DTM (dual wavelength approach). All images were averaged to 21 m pixel spacing as required by the Marrakech Accords to the Kyoto protocol. Color scale : 5…400 tC/ha.

236

H. Balzter et al. / Remote Sensing of Environment 108 (2007) 224–239

estimates is reduced to 36.0 tC/ha. A graph of all carbon content estimates for the 5 stands with 95% confidence intervals based on their standard deviation within each stand is given in Fig. 10. The figure shows that for the entire woodland the mean carbon content is similar between all methods, but that there is high between-stand variation in terms of the accuracy of the different methods. 3.3. Uncertainties in the carbon estimation Patenaude et al. (2003) described how a full quantification of the uncertainty in the allometric carbon accounting procedure cannot be achieved because of knowledge limitations. Therefore the impact of uncertainty in the allometric conversion factors on the canopy height to carbon conversion was assessed here using a Monte Carlo experiment. In the experiment all conversion factors were defined as Gaussian random variables where the standard deviations indicate the assumed uncertainties. The uncertainty propagation during the canopy height to carbon conversion was assessed for 1000 simulations, where the standard deviations in the Monte Carlo experiment were fixed at 10% of the value for the canopy height to top height conversion factor, biomass expansion factor, dry mass conversion factor, carbon content conversion factor, the constants for top height to stem volume conversion (Eq. (8)), and understorey carbon conversion factors (Eq. (9)). The random variation introduced to the mean carbon content estimate as a function of canopy height is shown in Fig. 11. The modeled carbon content increases exponentially as a function of canopy height as a result of the allometric equations used in the carbon accounting. As trees grow the root carbon mass and stem diameter increase in such a way that the overall carbon content is increasing non-linearly as a function of tree height. Younger trees grow in height more than in stem diameter, but older trees can grow in diameter more than in height. In contrast to pixel-based methods, clearly stand-

Fig. 10. Mean aboveground carbon content [t/ha] of 5 sampled forest stands at Monks Wood NNR based on the field survey (“TOTAL”), the LIDAR CHM, the XVV InSAR CHM using the XVV DSM and LIDAR DTM (“XVV”) and the dual-wavelength InSAR CHM (“XVVLHH”). The 95% confidence intervals show two standard deviations within each stand.

Fig. 11. Modeled aboveground carbon content [t/ha] determined from 1000 Monte Carlo simulations, where the allometric conversion factors were drawn from Gaussian distributions with defined uncertainty bounds. The 95% confidence interval shows twice the standard error of the mean carbon content from the 1000 runs.

based estimates reduce the standard error of the mean by averaging over an area. The achievable accuracy of stand-based carbon estimates increases with stand size. 4. Conclusions The penetration depth of microwaves into a forest canopy depends on the size and density distribution of the scattering elements in the resolution cell, the geometric structure of the scatterers and moisture conditions, as well as surface roughness and moisture content of the ground layer. Which components of a vegetation layer act as scattering elements is dependent on the wavelength and to a lesser extent the polarization characteristics of the SAR sensor. In the Monks Wood E-SAR dataset the scattering phase centers at L-band HH polarization were located at a range of different vertical heights ranging from the terrain underneath sparsely vegetated areas and in open areas such as rides and bridleways, to positions higher up in the canopy but below the canopy top. A spatial subset of high-coherence pixels was identified in which the scattering processes originated from a position sufficiently close to the terrain surface to allow the derivation of a smoothed interpolated DTM that includes both forested and unvegetated areas. Generally, in the L-HH DTM for dense forest areas the values are likely to stem from interpolations, while for sparse and open areas they tend to originate from the interferometric phase measurement. Further improvements in estimating the true terrain height from LHH are achievable, e.g. through using a forest structure model to account for the double-bounce contribution to the signal, which originates from trunk-ground interactions and may lead to an overestimation of the terrain height in some areas. Such an overestimation can be observed in Fig. 4, but it remains to be investigated whether double-bounce scattering is the cause for these observations. Martinez et al. (2000) found that forest vertical and horizontal heterogeneity must be addressed to fully understand X band backscatter mechanisms. Similar conclusions can be drawn about L-band from the study presented here.

H. Balzter et al. / Remote Sensing of Environment 108 (2007) 224–239

Limitations of the terrain height interpolation from L-band SAR interferometry are the dependence on vegetation canopy density, moisture conditions in the vegetation and ground layer and wave attenuation, the interpolation over areas with insufficient ground contribution to the signal and the phase noise. The DTM interpolation method uses a coherence threshold. If not enough pixels can be identified as yielding near-terrain height estimates, the algorithm smoothes out microscale topographic variability and leads to increasing errors. This is likely to be the case for more densely growing forests. There is a tradeoff between the coherence threshold in the interpolation step and the terrain model accuracy. With a high threshold a smaller number of high quality pixels are retained, leading to the loss of small scale topographic variation. With a lower threshold more pixels with poorer quality and ultimately contaminated with vegetation scattering are retained, also leading to a poorer spatial interpolation accuracy. The choice of an optimum coherence threshold is not trivial and is an issue for generalization of the method. This paper presents a case study for one test site, and the site-specific conditions are likely to limit the general applicability of the method. In some respects the test site Monks Wood is an unusual test sites for a radar remote sensing study in that it is a semi-natural deciduous woodland with a diverse and spatially heterogeneous vegetation canopy, which includes relatively sparse patches within the main wood. In contrast, for dense tropical rain forest Rignot (1996) derived canopy height from L-band InSAR where the observed scattering phase centers at L-band were closer to the canopy top than to the terrain. For dense vegetation types other ways of deriving a DTM have to be considered, e.g. from interferometric measurements over gaps in the canopy or from even longer wavelengths like P-band. A P-band SAR mission called BIOMASS is currently being studied by the European Space Agency as a possible future Earth Explorer mission. A polarimetric interferometric coherence optimization technique (Cloude & Papathanassiou, 1998) was developed to estimate terrain height under forest canopies. Its applicability under a range of different vegetation densities and structural forest types remains to be investigated. In summary, this study evaluates the retrieval accuracy of canopy height and derived carbon stock estimates from a new approach using dual-wavelength SAR interferometry from a combination of an X-band DSM and a smoothed interpolated DTM derived from L-HH interferometry over high-coherence pixels. The error figures (rmse) of carbon content per hectare contain error components from the canopy height estimation (relative error = 28.5% when compared to the LIDAR canopy height model, which has itself a relative error of 17.6% compared to theodolite data), the allometric conversion to carbon content (unknown uncertainty), and from the independent field survey that provided a sample-based estimate of carbon stocks per forest stand (relative error is thought to exceed 10% but cannot be fully quantified, see Patenaude et al., 2003 for details). The dual-wavelength InSAR technique could be applied to data from a constellation like the TerraSAR-L and Tandem-X missions, whereas imaging LIDAR systems are likely to rely on airborne platforms in the near future. Sensors

237

like ICESAT-GLAS (Lefsky et al., 2005) or the proposed Vegetation Canopy LIDAR (VCL) provide a much coarser spatial sampling with good data on vertical profiles of vegetation canopies, but higher uncertainties about the spatial distribution of vegetation. The necessary adjustment of the dualwavelength InSAR CHM for penetration into the canopy through gaps in the foliage is likely to be test site specific, and is a limitation for a generalization of the method, unless an algorithm to model or estimate the penetration depth from the observations can be found. The rmse of the dual-wavelength InSAR CHM could also potentially be reduced if single-pass acquisitions were used for both wavelengths. A comprehensive coherent microwave modeling experiment is needed to quantify the expected variability of the L-HH scattering phase center height as a function of vegetation type and structure, soil and canopy moisture content and temporal acquisition lag. It should include an assessment of the sensitivity of X-band interferometry to variations in canopy structure and moisture. A single-pass interferometric acquisition at X-band is critical to minimize phase noise. Potential future observing systems like the single-pass interferometric system TANDEM-X that would add value to TerraSAR-X would be able to deliver ~ 1 m spatial resolution interferometric elevation data of the Earth's surface at an accuracy of about 2 m. In combination with a repeat-pass L-band interferometer based on TerraSAR-L or ALOS it may be possible to evaluate the technique over a set of global test sites to quantify the reliability and limitations of the method, and particularly to delimit the range of densities of scattering elements in different forest ecosystems with varying canopy structure conditions for which the algorithm gives reliable terrain estimates. Longer wavelengths such as P-band generally penetrate deeper into dense forest canopies than Lband, and are likely to yield a larger number of pixels where the radar return originates from near the terrain surface under such vegetation. Acknowledgements This study was supported by the UK Natural Environment Research Council (NERC) under the New Observing Techniques Programme, CORSAR project, Grant number NER/Z/S/ 2000/01282. ©E-SAR data by British National Space Center and NERC. LIDAR data are courtesy of the Environment Agency, our thanks go to Neil Veitch for providing the data. The work was carried out at the Center for Ecology and Hydrology Monks Wood, UK. Ross Hill and David Gaveau (Center for Ecology and Hydrology, Monks Wood) processed and validated the LIDAR data. Urs Wegmüller (Gamma Remote Sensing AG, Bern) provided valuable discussions about interferometric data processing. Terry Dawson (University of Edinburgh), Laine Skinner and Adrian Luckman (University of Swansea) supported the field data collection. Genevieve Patenaude (University of Edinburgh) developed the carbon accounting equations. Rody Nigel (University of Cambridge) provided support with the equations in the methods section. We thank the reviewers for their constructive comments which have helped to improve the manuscript.

238

H. Balzter et al. / Remote Sensing of Environment 108 (2007) 224–239

References Andersen, H. -E., McGaughey, R., Reutebuch, S., Schreuder, G., Agee, J., & Mercer, B. (2004). Estimating canopy fuel parameters in a Pacific Northwest conifer forest using multifrequency polarimetric IFSAR. Proceedings of the ISPRS conference, Istanbul. Askne, J. I. H., Dammert, P. B. G., Ulander, L. M. H., & Smith, G. (1997). Cband repeat-pass interferometric SAR observations of the forest. IEEE Transactions on Geoscience and Remote Sensing, 35, 25−35. Askne, J., & Santoro, M. (2005). Multitemporal repeat pass SAR interferometry of boreal forests. IEEE Transactions on Geoscience and Remote Sensing, 43, 1219−1228. Aziz, H. K., & White, K. (2003). Using mimics to model L-band SAR backscatter from a peat swamp forest. Journal of Tropical Forest Science, 15, 546−556. Balzter, H. (2001). Forest mapping and monitoring with interferometric Synthetic Aperture Radar (InSAR). Progress in Physical Geography, 25, 159−177. Balzter, H., Baker, J. R., Hallikainen, M., & Tomppo, E. (2002). Retrieval of timber volume and snow water equivalent over a Finnish boreal forest from airborne polarimetric Synthetic Aperture Radar. International Journal of Remote Sensing, 23, 3185−3208. Balzter, H., Saich, P., Luckman, A. J., Skinner, L., & Grant, J. (2002). Forest stand structure from airborne polarimetric InSAR. In Ag Noordwijk (Ed.), Proceedings of the Third International Symposium on retrieval of bio- and geophysical parameters from SAR data for land applications (pp. 321−326). ESA Publications Division C/O ESTEC. Balzter, H., Skinner, L., Luckman, A., & Brooke, R. (2003). Estimation of tree growth in a conifer plantation over 19 years from multi-satellite L-band SAR. Remote Sensing of Environment, 84, 184−191. Bamler, R., Adam, N., Davidson, G. W., & Just, D. (1998). Noise-induced slope distortion in 2-D phase unwrapping by linear estimators with application to SAR interferometry. IEEE Transactions on Geoscience and Remote Sensing, 36, 913−921. Bamler, R., & Hartl, P. (1998). Synthetic Aperture Radar interferometry. Inverse Problems, 14, R1−R54. Cloude, S. R., & Papathanassiou, K. P. (1998). Polarimetric SAR interferometry. IEEE Transactions on Geoscience and Remote Sensing, 36, 1551−1565. Costantini, M. (1998). A novel phase unwrapping method based on network programming. IEEE Transactions on Geoscience and Remote Sensing, 36, 813−821. Costantini, M., Farina, A., & Zirilli, F. (1999). A fast phase unwrapping algorithm for SAR interferometry. IEEE Transactions on Geoscience and Remote Sensing, 37, 452−460. Eriksson, L. E. B., Santoro, M., Wiesmann, A., & Schmullius, C. C. (2003). Multitemporal JERS repeat-pass coherence for growing-stock volume estimation of Siberian forest. IEEE Transactions on Geoscience and Remote Sensing, 41, 1561−1570. Forestry Commission (1981). Yield models for forest management (booklet 48). Edinburgh: Forestry Commission. Fornaro, G., & Sansosti, E. (1999). A two-dimensional region growing least squares phase unwrapping algorithm for interferometric SAR processing. IEEE Transactions on Geoscience and Remote Sensing, 37, 2215−2226. Fornaro, G., Franceschetti, G., & Lanari, R. (1996). Interferometric SAR phase unwrapping using Green's formulation. IEEE Transactions on Geoscience and Remote Sensing, 34, 720−727. Fornaro, G., Franceschetti, G., Lanari, R., Rossi, D., & Tesauro, M. (1997). Interferometric SAR phase unwrapping using the finite element method. IEE Proceedings. Radar, Sonar and Navigation, 144, 266−274. Gaveau, D. (2002). Modelling the dynamics of ERS-1/2 coherence with increasing woody biomass over boreal forests. International Journal of Remote Sensing, 23, 3879−3885. Griffiths, H. D., & Wilkinson, A. J. (1994). Improvements in Phase Unwrapping Algorithms for Interferometric SAR. Onde Electrique, 74, 46−52. Hagberg, J. O., Ulander, L. M. H., & Askne, J. (1995). Repeat-Pass SAR Interferometry Over Forested Terrain. IEEE Transactions on Geoscience and Remote Sensing, 33, 331−340. Hamilton, G. J. (1975) Forest Mensuration. Forestry Commission Booklet no. 39. The Stationary Office, London, HMSO, 274 pp.

Hill, R. A., Gaveau, D. L. A., & Spendlove, M., 2002, Accuracy issues in the assessment of deciduous woodland canopy height by airborne laser scanning: a case study. Proc. of ForestSAT, 5-9 August 2002, Edinburgh, Forest Research, Forestry Commission. Hill, R. A., & Thomson, A. G. (2005). Mapping woodland species composition and structure using airborne spectral and LiDAR data. International Journal of Remote Sensing, 26, 3763−3779. Houghton, R. A. (2005). Aboveground forest biomass and the global carbon balance. Global Change Biology, 11, 945−958. Hyyppä, J., & Hallikainen, M. (1996). Applicability of airborne profiling radar to forest inventory. Remote Sensing Environment, 57, 39−57. Imhoff, M. L. (1995). Radar backscatter and biomass saturation- ramifications for global biomass inventory. IEEE Transactions on Geoscience and Remote Sensing, 33, 511−518. Kasischke, E. S., Melack, J. M., & Dobson, M. C. (1997). The use of imaging radars for ecological applications — a review. Remote Sensing of Environment, 59, 141−156. Kellndorfer, J., Walker, W., Pierce, L., Dobson, C., Fites, J. A., Hunsaker, C., et al. (2004). Vegetation height estimation from shuttle radar topography mission and national elevation datasets. Remote Sensing of Environment, 93, 339−358. Le Toan, T., Beaudoin, A., Riom, J., & Guyon, D. (1992). Relating Forest Biomass to SAR Data. IEEE Transactions on Geoscience and Remote Sensing, 30, 403−411. Lefsky, M. A., Harding, D. J., Keller, M., Cohen, W. B., Carabajal, C. C., Espirito-Santo, F. D., et al. (2005). Estimates of forest canopy height and aboveground biomass using ICESat. Geophysical Research, 32. Lin, Q., Vesecky, J. F., & Zebker, H. A. (1994). Phase unwrapping through fringe-line detection in synthetic-aperture radar interferometry. Applied Optics, 33, 201−208. Lyuboshenko, I., & Maitre, H. (1999). Phase unwrapping for interferometric Synthetic Aperture Radar by use of Helmholtz equation eigenfunctions and the first Green's identity. Journal of the Optical Society of America. A, Optics, Image Science, and Vision, 16, 378−395. Martinez, J. M., Floury, N., Le Toan, T., Beaudoin, A., Hallikainen, M., & Makynen, M. (2000). Measurements and modeling of vertical backscatter distribution in forest canopy. IEEE Transactions on Geoscience and Remote Sensing, 2000, 710−719. Madsen, S. N., & Zebker, H. A. (1998). Imaging Radar Interferometry. Principles and Applications of Imaging Radar, Manual of Remote Sensing, 3rd Edition Soc. of Photogrammetry and Rem. Sens., Vol. 2. (pp. 359−406) Chapter 6. Papathanassiou, K. P., & Cloude, S. R. (2001). Single-baseline polarimetric SAR interferometry. IEEE Transactions on Geoscience and Remote Sensing, 39, 2352−2363. Patenaude, G. L., Briggs, B. D. J., Milne, R., Rowland, C. S., Dawson, T. P., & Pryor, S. N. (2003). The carbon pool in a British semi-natural woodland. Forestry, 76, 109−119. Patenaude, G. L., Hill, R. A., Milne, R., Rowland, C. S., & Dawson, T. P. (2002). Forest carbon accounting using Airborne Laser Scanning remote sensing and modelling approaches. Proc. of ForestSAT 2002, 5-9 August Edinburgh, UK: Forest Research. Patenaude, G., Hill, R. A., Milne, R., Gaveau, D. L. A., Briggs, B. B. J., & Dawson, T. P. (2004). Quantifying forest above ground carbon content using LiDAR remote sensing. Remote Sensing of Environment, 93, 368−380. Patenaude, G., Milne, R., & Dawson, T. P. (2005). Synthesis of remote sensing approaches for forest carbon estimation: reporting to the Kyoto Protocol. Environmental Science & Policy, 8, 161−178. Pritt, M. D. (1996). Phase unwrapping by means of multigrid techniques for interferometric SAR. IEEE Transactions on Geoscience and Remote Sensing, 34, 728−738. Prusinkiewicz, P., & Lindenmayer, A. (1990). The algorithmic beauty of plants. New York: Springer. Pulliainen, J., Engdahl, M., & Hallikainen, M. (2003). Feasibility of multitemporal interferometric SAR data for stand-level estimation of boreal forest stem volume. Remote Sensing of Environment, 85, 397−409. Rignot, E. (1996). Dual-frequency interferometric SAR observations of a tropical rain-forest. Geophysical Research Letters, 23, 993−996.

H. Balzter et al. / Remote Sensing of Environment 108 (2007) 224–239 Rignot, E., Way, J. B., Williams, C., & Viereck, L. (1994). Radar estimates of aboveground biomass in boreal forests of interior alaska. IEEE Transactions on Geoscience and Remote Sensing, 32, 1117−1124. Rosenqvist, A., Milne, A., Lucas, R., Imhoff, M., & Dobson, C. (2003). A review of remote sensing technology in support of the Kyoto Protocol. Environmental Science & Policy, 6, 441−455. Saich, P., Balzter, H., Luckman, A., & Skinner, L. (2001). Deriving forest characteristics using polarimetric InSAR measurements and models. Proceedings of the International Geoscience and Remote Sensing Symposium, IGARSS 2001, Sydney, 9-13 July CDROM. Santoro, M., Askne, J., & Dammert, P. B. G. (2003). Tree height estimation from multi-temporal ers sar interferometric phase. Proc. Fringe, ESA-ESRIN, 1–5 Dec. Schimel, D. S., House, J. I., Hibbard, K. A., Bousquet, P., Ciais, P., Peylin, P., et al. (2001). Recent patterns and mechanisms of carbon exchange by terrestrial ecosystems. Nature, 414, 169−172. Stebler, O., Meier, E., & Nuesch, D. (2002). Multi-baseline polarimetric SAR interferometry — first experimental spaceborne and airborne results. ISPRS Journal of Photogrammetry and Remote Sensing, 56, 149−166. Tansey, K. J., Luckman, A. J., Skinner, L., Balzter, H., Strozzi, T., & Wagner, W. (2004). Classification of forest volume resources using ERS tandem coherence and JERS backscatter data. International Journal of Remote Sensing, 25, 751−768. Treuhaft, R. N., Asner, G. P., & Law, B. E. (2003). Structure-based forest biomass from fusion of radar and hyperspectral observations. Geophysical Research, 30.

239

Treuhaft, R. N., Law, B. E., & Asner, G. P. (2004). Forest attributes from radar interferometric structure and its fusion with optical remote sensing. Bioscience, 54, 561−571. Treuhaft, R. N., Madsen, S. N., Moghaddam, M., & vanZyl, J. J. (1996). Vegetation characteristics and underlying topography from interferometric radar. Radio Science, 31, 1449−1485. Trouve, E., Nicolas, J. M., & Maitre, H. (1998). Improving phase unwrapping techniques by the use of local frequency estimates. IEEE Transactions on Geoscience and Remote Sensing, 36, 1963−1972. Wang, Z. J., & Li, S. S. (1999). Phase unwrapping through a branch-cut-based cut-bridging and window-patching method. Applied Optics, 38, 805−814. Wegmüller, U., & Werner, C. (1997). Retrieval of vegetation parameters with SAR interferometry. IEEE Transactions on Geoscience and Remote Sensing, 35, 18−24. Xu, W., & Cumming, I. (1999). A region-growing algorithm for InSAR phase unwrapping. IEEE Transactions on Geoscience and Remote Sensing, 37, 124−134. Zebker, H. A., & Lu, Y. P. (1998). Phase unwrapping algorithms for radar interferometry: Residue-cut, least-squares, and synthesis algorithms. Journal of the Optical Society of America. A, Optics, Image Science, and Vision, 15, 586−598. Zebker, H. A., & Villasenor, J. (1992). Decorrelation in Interferometric Radar Echoes. IEEE Transactions on Geoscience and Remote Sensing, 30, 950−959.