Tropical-Forest Structure and Biomass Dynamics from TanDEM-X Radar Interferometry

Changes in tropical-forest structure and aboveground biomass (AGB) contribute directly to atmospheric changes in CO2, which, in turn, bear on global climate. This paper demonstrates the capability of radar-interferometric phase-height time series at X-band (wavelength = 3 cm) to monitor changes in vertical structure and AGB, with sub-hectare and monthly spatial and temporal resolution, respectively. The phase-height observation is described, with a focus on how it is related to vegetation-density, radar-power vertical profiles, and mean canopy heights, which are, in turn, related to AGB. The study site covers 18× 60 km in the Tapajós National Forest in the Brazilian Amazon. Phase-heights over Tapajós were measured by DLR’s TanDEM-X radar interferometer 32 times in a 3.2 year period from 2011–2014. Fieldwork was done on 78 secondary and primary forest plots. In the absence of disturbance, rates of change of phase-height for the 78 plots were estimated by fitting the phase-heights to time with a linear model. Phase-height time series for the disturbed plots were fit to the logistic function to track jumps in phase-height. The epochs of clearing for the disturbed plots were identified with ≈1-month accuracy. The size of the phase-height change due to disturbance was estimated with ≈2-m accuracy. The monthly time resolution will facilitate REDD+ monitoring. Phase-height rates of change were shown to correlate with LiDAR RH90 height rates taken over a subset of the TanDEM-X data’s time span (2012–2013). The average rate of change of phase-height across all 78 plots was 0.5 m-yr−1 with a standard deviation of 0.6 m-yr−1. For 42 secondary forest plots, the average rate of change of phase-height was 0.8 m-yr−1 with a standard deviation of 0.6 m-yr−1. For 36 primary forest plots, the average phase-height rate was 0.1 m-yr−1 with a standard deviation of 0.5 m-yr−1. A method for converting phase-height rates to AGB-rates of change was developed using previously measured phase-heights and field-estimated AGB. For all 78 plots, the average AGB-rate was 1.7 Mg-ha−1-yr−1 with a standard deviation of 4.0 Mg-ha−1-yr−1. The secondary-plot average AGB-rate was 2.1 Mg-ha−1-yr−1, with a standard deviation of 2.4 Mg-ha−1-yr−1. For primary plots, the AGB average rate was 1.1 Mg-ha−1-yr−1 with a standard deviation of 5.2 Mg-ha−1-yr−1. Given the standard deviations and the number of plots in each category, rates in secondary forests and all forests were significantly different from zero; rates in primary forests were consistent with zero. AGB-rates were compared to change models for Tapajós and to LiDAR-based change measurements in other tropical forests. Strategies for improving AGB dynamical monitoring with X-band interferometry are discussed. Forests 2017, 8, 277; doi:10.3390/f8080277 www.mdpi.com/journal/forests Forests 2017, 8, 277 2 of 28


Introduction
Tropical forests, which account for about 50% of the world's forested biomass, play a critical role in the control of atmospheric carbon dioxide through emissions from forest removals and uptake through forest accruals [1].Some hypotheses suggest that some of the missing carbon sink is due to uptake in mature tropical forests [2,3].Monitoring of changes in aboveground biomass (AGB) in tropical forests at the fine spatial scales at which disturbance takes place, 25-100 m [4], is required for global carbon cycle characterization.Climate change mitigation efforts designed to reduce emissions from deforestation and forest degradation (REDD+) will require monitoring much more frequently than the ≈1 time per year realized by optical satellites in the tropics [5].Monitoring on the temporal scales at which disturbance takes place, ≈1 month, will be required to enable REDD+ [6].
Interferometric SAR (InSAR), with its broad, all-weather coverage, and sensitivity to vegetation vertical structure [7,8], has the potential to play a key role in AGB-change monitoring strategies at local to regional to global scales.We present interferometric phase-height change from the 2-spacecraft, X-band TanDEM-X system [9], over Tapajós National Forest in Brazil.Interferometric phase-height is related to vertical-profile-averaged height, as described in Section 2.3.Phase-heights were measured up to 15 times per year from 2011-2014 on 0.25-ha plots.The high temporal density of measurements was possible because, unlike optical signals, radar signals at X-band are not attenuated by clouds and rain.Unlike radar power measurements, InSAR measurements of phase-height and coherence are directly sensitive to vertical structure, from which AGB can be estimated.(See discussion of (1) in Section 2.2 for definitions of "phase-height" and "coherence").While X-band power has proven inaccurate for forest AGB measurement, InSAR phase-height and InSAR coherence produce reasonable (≈30%, 57 Mg-ha −1 ) single-epoch AGB estimates for tropical forests at X-band [10].At Tapajós the rates for various LiDAR metrics were available over the 2012 and 2013 period and were compared to the InSAR phase-height rates.Using time series of phase-heights, we found that we could estimate changes in AGB 2-4 times more accurately than we could estimate single-epoch, static AGB from a single measurement of phase-height.
The purpose of this paper is to assess and demonstrate the potential accuracy of rates of change of tropical-forest structure and AGB estimation available from an X-band interferometer with aggregated spatial resolution of 0.25 ha and temporal resolution of ≈1 month.One previous study used 10 observation epochs of TanDEM-X data spanning 2 years in the tropics to gauge the stability of an 800 m × 800 m area [11], but the referenced study showed stability levels of the large area, rather than rates of change at the ha scale.Other studies [12,13] used multi-temporal TanDEM-X data in boreal forests largely to demonstrate stability rather than change.Another study used X-band InSAR from SRTM and TanDEM-X to measure forest change over an 11-year period, 2000 to 2011 [14], in boreal forests.Studies have been done which use TanDEM-X InSAR data to estimate AGB in the tropics at single-epochs [10,15,16].One study combined LiDAR and TanDEM-X interferometry [17], but, again for single-epoch, and for temperate forests.LiDAR has been used to infer AGB change from chronosequences in the tropics [18].AGB tropical dynamics measured directly with LiDAR have been addressed by [19,20].However, as opposed to this study which had multiple observation epochs per year, these two LiDAR studies had two observations each, separated by 7 and 10 years in La Selva, Costa Rica, and Barro Colorado Island, Panama, respectively.Thus the current study using TanDEM-X phase-height seems to fill a gap of reporting significant structure and AGB-rates in both primary and secondary tropical forests with a temporal resolution of ≈1 month.
An ancillary objective of this paper is to describe a general framework for transforming phase-height rate to AGB-rate.Various semi-empirical expressions exist to give standing, single-epoch AGB in terms of height, Lorey's height, or phase-height, e.g., [12,16,21], but we describe general Forests 2017, 8, 277 3 of 28 considerations for transforming phase-height rates to AGB-rates.While we use a power law to describe the mass density function in (11), the formalism described shows how any functional form for the density function could be accommodated.
Section 2 describes Tapajós National Forest and TanDEM-X.It relates the InSAR phase-height to vertical characteristics of the forest, and to the power-averaged mean canopy height, which is defined in (2) in Section 2.3.After discussing the estimation and calibration of phase-height data from TanDEM-X raw data, Section 3 shows phase-height time series results over the 2011-2014 period for secondary and primary forest plots.It further shows rates of phase-heights from linear fits over time, as a function of field-estimated AGB, suggesting that TanDEM-X phase-height can distinguish between the faster growth of secondary plots, and the slower rate of change of old growth.Phase-height rates over one year-August 2012-August 2013-are compared to the rates of LiDAR's RH90-the height below which 90% of the LiDAR power is returned, starting from the peak of the ground signal.A conversion factor is derived to convert phase-height rates of change to AGB-rates of change, and AGB-rates are calculated for the 78 plots.Section 4 is a discussion of the results, the choice of phase-height rate to estimate AGB rate, error analysis, and future improvements.

Materials and Methods
This section describes Tapajós National Forest, field vegetation measurements, the InSAR phase-height, and how it is related to radar-power profiles, and the power-averaged mean in particular.It also describes the LiDAR data used for verification of InSAR phase-height rate.

Tapajós National Forest and Field Measurements
The Tapajós National Forest is ≈50 km south of Santarem, Pará, in the central Brazilian Amazon.Its climate, rainfall, and vegetation are as described in [10].Field measurements were taken in September 2010 on 26 plots (50 m × 50 m = 0.25 ha) and on 52 plots of the same size in August 2013.The maximum value of AGB was 738 Mg-ha −1 , and averaged 189 Mg-ha −1 .The plots were primary and secondary forests, and included some selective logging sites.Througout this paper, "primary" means a forest that has never been logged and has developed following natural disturbances and under natural processes."Secondary" means a forest developing after a stand-replacing disturbance (i.e., a disturbance that eliminated all, or most, previous trees in the stand) such as logging or the clearing of land for agriculture or cattle.Figure 1 shows the locations of 78 plots along with the large rectangle corresponding to TanDEM-X coverage used in this study, and the small red rectangles indicating LiDAR coverage in 2012 and 2013.The field plots were representative of the variety of the area observed by TanDEM-X and LiDAR.Field measurements for each tree included total height, height-to-base-of-crown, tree diameter at breast height, species (to assign wood density values derived from the literature), crown dimensions, and location.The crown dimensions were the distance between extrema of leaf area as viewed from the ground, along two orthogonal axes [22].The diameter, height, and wood density were used with allometric equations [23] to estimate AGB.The accuracy of the AGB estimates was 25% [24].Plots were geolocated with sub-meter accuracy using differential GPS and a total station.

TanDEM-X, The Interferometric Phase-Height and Its Relation to Radar Power Vertical Profiles
Built and operated by the German Aerospace Center (DLR), the TanDEM-X interferometer is the only 2-spacecraft interferometer to orbit the Earth to date [9].Depending on the technical objective (e.g., forest research, digital elevation map), the vector separation of the two spacecraft, called the "baseline", can be adjusted to have lengths under 100 m up to several kilometers.For our forest data over Tapajós, baseline separations ranged from 50 m to 300 m, but were mostly in the 180-m range.All TanDEM-X data in this paper were acquired at horizontal polarization, both transmit and receive.The incidence angle was approximately 40 • .
The phase-height, the principal indicator of change in this paper, results from the interferometric phase of forests, which have vertically distributed vegetation.In order to understand the phase of a forest with many vertically distributed scatterers, first consider the interferometric measurement of a single scatterer, at a single altitude, say, the top of the tree in Figure 2. The ends of the TanDEM-X interferometer labeled 1 and 2 in Figure 2 receive the signals E 1 and E 2 , which propagate along the two paths, of lengths r 1 and r 2 .The single-scatterer interferometric phase is proportional to the difference r 1 -r 2 , which, in turn, depends on the height of the scatterer; thus the InSAR phase depends on the height of the scatterer.The single-scatterer interferometric phase is obtained by cross multiplying E 1 and E * 2 , where the star is the complex conjugate.For forests consisting of many vertically distributed scatterers (leaves, branches, etc.) the interferometric phase results from all scatterers, and is obtained from the complex cross correlation of E 1 and E * 2 as in (1) below.In order to show the relationship between interferometric phase-height and the radar power profile p(z ) for the forest, the cross correlation of the E 1 and E 2 signals is expressed as proportional to the sum (integral) of the contributions from scatterers at z', extending from the ground up to a forest height of h v [7]: where the ensemble average brackets indicate spatial averaging over scatterer positions and strengths from different small (few-meter) looks to aggregate to the plot size-0.25 ha in this study.The integral can be seen as summing the contributions at each altitude z of a phasor (complex exponential) to the InSAR cross correlation with phase α z z .The constant α z is the derivative of InSAR phase with respect to altitude, also known as the "vertical wavenumber".It is proportional to the length of the baseline, B, and other geometric factors.It can also be thought of as the Fourier frequency over which < p(z ) > is Fourier transformed in (1).The integral's aggregate amplitude A and phase, φ int , are indicated.The phase-height, h φ , is the single-scatterer height which would produce the same phase as that of the entire forest, and, as indicated in (1), h φ = φ int /α z .The coherence is A divided by the same integral with α z = 0.The often-used quantity, the ambiguity height, h amb ≡ 2π/α z , is the height at which interferometric phase winds through a cycle.For many of the measurements in this report, h amb ≈ 80 m, meaning the top and bottom of an 80 m tree would have the same phase (0 • = 360 • ).Trees of the order of the ambiguity height would complicate interpretations of phase, so, for the most part, only h amb taller than the tallest trees in Tapajós (≈45 m), which are of utmost importance in the structural arrangement of the forest, were used.A list of dates of passes with ambiguity heights is in Table A1.The amplitude of each phasor's contribution in (1) at z -the average radar power < p(z ) >-depends on three fundamental descriptors common to virtually all remote sensing: (1) the spatial density of scatterers-leaves, branches, trunks, ground, as a function of height, (2) the electromagnetic strengths of those scatterers, and (3) the attenuation of incident or scattered fields in the medium.

The Relationship of Phase-Height to Mean Canopy Height and AGB
Each of the three descriptors on which radar power depends has biological analogues.The spatial density of scatterers (number of scatterers per unit volume) is related to leaf-area-density (LAD), in that LAD depends on both the area of leaves and the number of leaves per unit volume.Qualitatively, the larger the values of the spatial density profile at high altitudes, contributing to p(z ), the taller the vegetation, and the higher the AGB, where the ensemble average brackets are understood from here on.The second descriptor contributing to p(z ), the brightness profile, bears on leaf area density also, as larger leaves have larger area and are frequently better backscatterers of the radar.The third descriptor, the attenuation, like the scatterer density, has a simple biological analog.Attenuation increases when scatterer/vegetation density increases, and/or when the water content of the leaves increases, both of which can conceivably correlate with AGB increase.Because radar-power-profile-averaged canopy height (defined in (2) below)-called "mean canopy height" or "MCH"-depends on the vertical distribution of scatterer density and brightness in p(z ), and because the distribution of scatterer density is related to AGB, it is plausible that AGB can be associated with some function of radar MCH.The MCH is and, in part for its relation to AGB stated above, we postulate that some function of the radar MCH is a good indicator of AGB.Another reason to postulate the utility of the radar MCH is based on the LiDAR MCH.It is given by (2), replacing p(z ) with the LiDAR waveform.Linear, logarithmic, or power-law functions of the LiDAR MCH have been used as significant AGB independent variables [20,[25][26][27].The success of the above-mentioned functions of LiDAR MCH to predict AGB in part prompts exploring functions of the radar MCH to monitor AGB change in this paper.If we therefore want to model changes in AGB based on changes in radar MCH but, from (1) we only have phase-height at our disposal, a natural question is, for a given p(z ), how close is phase-height from (1) to MCH from (2)?And how close are phase-height rates to MCH rates?Section 2.3.1 presents a qualitative argument that phase-height is close to MCH.Section 2.3.2 presents quantitative scenarios where MCH is exactly equal to phase-height, and Section 2.3.3 quantitatively investigates the general magnitude of the differences between phase-height and MCH, as well as differences in their rates.

Qualitative Relationship Between Phase-Height and Mean Canopy Height
In order to qualitatively address the proximity of phase-height to MCH, Figure 3 is a graphic representation of the integral in (1), for a profile p(z ) with just 3 thin vegetation layers at z 1 , z 2 , and z 3 .Then the integral in (1) becomes a sum of 3 vectors in the complex plane, with the length of each black vector given by the power coming from that layer, p(z 1 ), p(z 2 ), and p(z 3 ); and the phase of each layer-the angle between the vector and the x-axis-proportional to each layer's vertical position, α z z 1 , α z z 2 , and α z z 3 .The vector in red, with total phase, φ int , and amplitude, A, is the value of the integral.The total phase-height, φ int /α z = h φ , is between the lowest (z 1 ) and highest (z 3 ) values of the heights in the profile, while MCH is also between the extrema of the profile.Furthermore, the phasor corresponding to each phase-height z is vectorially added with a strength proportional to p(z ), while in the expression for MCH, (2), the term for each height z is linearly added with a strength proportional to p(z ). Figure 3 and the form of (1) and (2) therefore suggest the plausibility that phase-height is related to radar MCH.

Real
Imaginary A graphic representation of the integral in (1) for a 3-thin-layer radar power profile.
The ntegral is a sum of the three black vectors in the complex plane, each with amplitude equal to the power of the layer, and phase proportional to the altitude of the layer.The red vector is the total, the value of the integral in (1).The phase-height corresponding to this red vector, φ int /α z , is between that of the lowest and the highest layers.

Quantitative Relationship Between Phase-Height and Mean Canopy Height: Conditions for Equality
A more quantitative examination of the relationship between phase-height (h φ ) and MCH starts with the expression for MCH above, and, from (1), the phase-height: There are two conditions when h φ is exactly equal to MCH.The first is when α z h v is very small and only the first two terms of the Taylor expansion of the complex exponential are kept in (3): In order for α z to be very small, h amb must be very high.For a given level of phase noise, very high h amb (low α z ) will induce high levels of noise in phase-height, which is in part why using very large h amb is not always practical.
The second condition is if p(z ) is symmetric about some value of z , say z 0 : where g is any real function, as long as it depends only on the magnitude of the difference between z and the symmetry point z 0 .Inserting (5) into (2) confirms that MCH = z 0 , if the lower and upper limits of integration are extended to ±∞ and p(z ) constrains the integral.Inserting (5) into (1), extending the lower and upper limits to ±∞ again, the integral becomes proportional to which is a real integral multiplying a complex phasor.Only the complex phasor determines the phase, and therefore the phase-height, in (3), which is h φ = z 0 .The phase-height equals the symmetry height, which equals MCH.In general, the two conditions in ( 4) and ( 5) are not met, α z h v is not <1, and p(z ) is not perfectly symmetric; and the phase-height is not equal to MCH, nor are their rates equal.In order to get an estimate of the order of magnitude of the differences between h φ and MCH due to asymmetry in p(z ), we inserted the following model radar power profile into (2) and ( 3), which were numerically integrated to calculate MCH and h φ : where σ le f t and σ right are the left-and right-handed standard deviations.σ le f t was held at 2 m, and z 0 , the symmetry point when σ le f t = σ right , was set equal to 12.5 m for this calculation.A1).Note that for σ right = 2 m, on the far left of Figure 4, the difference h φ -MCH = 0 on the y axis of Figure 4, as in ( 6), for the symmetric profile.When asymmetry is introduced into p(z ) by letting σ right increase, the h φ − MCH difference departs from zero, and for the most asymmetric profile, that difference is ≈−0.6 m for h amb = 60 m and ≈−0.035 m for h amb = 250 m.Because h amb 's for different epochs were provided by the TanDEM-X mission, the differences in h φ in a time series arise in part due to different departures of h φ from MCH owing to h amb diversity.It is therefore the difference in the curves in Figure 4 that constitutes an estimate of the error in using h φ instead of MCH.That is, if over the 3.2-year observation period, forest structure of a particular plot did not change, 0.6 m is a worst-case estimate of the RMS scatter of the fit of phase-height to time, due to changes in observation baseline (and h amb ) causing changes in h φ .This scatter would be due to changes in the observing instrument and not in the forest itself.Similarly, from Figure 5 which gives phase-height rate errors and RMS scatters about a linear fit to time, a 0.6 m scatter for all points would generate about 0.1 m-yr −1 rate error in using h φ instead of MCH.These are worst-case estimates for volume asymmetries because the actual temporal distribution of ambiguity heights did not begin and end with the extremes of 50 m and 250 m.Rather these extremes were imbedded in a more even 70-140 m range (Table A1).h amb 's will therefore not generally differ from each other as much as the extreme case used for illustration in Figure 4. Also, the most extreme asymmetric profile (right-hand inset of p(z') versus altitude above the ground) leads to the magnitude 0.6 m h φ -MCH difference.Actual p(z )s will probably not be as asymmetric as the right-hand inset, though this should be checked with radar and LiDAR profiles [22].That 0.6 m and 0.1 m-yr −1 are smaller than the scatter and rate error observed in fits of phase-height to time, and that these values are worst case, suggests that phase-height can be taken as nearly the same as MCH for phase-height rate Forests 2017, 8, 277 9 of 28 analysis.Excluding very high h amb from the analysis did not change results very much, and this is another suggestion that the order of magnitude of the differences between phase-height and MCH is at the meter level.However, the simple test with (7) does not treat possible asymmetric ground contributions, under the assumption that at X-band, we are mostly seeing the volume.The key point of Figure 4 is that asymmetry in p(z ) is the major contributor to differences between phase-height and MCH, and volume asymmetry causes a meter-level order-of-magnitude difference between MCH and h φ , with ≈0.1 m-yr −1 differences in rates.Note that in a future mission with forest structure and biomass dynamics as a priority, this source of error in using h φ could be greatly reduced by using the same h amb for all observation epochs.Then, for the h amb used, only the derivative of the curve like those in Figure 4, dMCH is relevant.The derivative gives the change in the difference between h φ and MCH due to changes in MCH, i.e. due to changes in forest structure.Analysis of this derivative for h amb = 60 m suggests that less than ≈0.1 m errors result from using h φ relative to using the radar MCH.The calculated difference between phase-height and mean canopy height (power-profile-averaged height), realized for an asymmetric-Gaussian hypothetical model of the volume radar power profile p(z ) as in (7).σ le f t is held at 2 m.When σ le f t = σ right , p(z ) is symmetric, the h φ -MCH difference is zero, and MCH = 12.5 = z 0 .The symmetric profile of p(z ) versus altitude in meters is in the left-hand inset.With σ right = 12, the right-hand inset shows the highly asymmetric power profile as a function of altitude in meters, and the h φ -MCH difference is −0.6 m.The h φ -MCH difference is plotted against MCH, which goes from 12.5 to 20 m when σ right is increased from 2 to 12 m.The difference is shown for two heights of ambiguity.

LiDAR Data for Phase-Height Rate Verification
In the next section, rates of LiDAR metrics were compared to phase-height rates.LiDAR data available for 25 of the plots were taken by GEOID Laser Mapping Ltda. in 2012 and 2013.They were acquired in 2012 using an ALTM 3100 (Teledyne Optech, 300 Interchange Way Vaughan, ON, Canada), and data acquired in 2013 used ALTM Orion M-200 (Teledyne Optech).Instruments were flown at 850-900 m above the ground, with swath side overlap of 65% and a field of view of 11 • .For each LiDAR pulse, up to 4 returns were recorded, with approximately 25-39 returns m −2 .A minimum return density of ≥4 m −2 was realized over 99.5% of the LiDAR area.

Results: Time Series of Phase-Height and Aboveground Biomass
This section describes how phase-heights are derived from raw InSAR data.It also describes how systematic trends not related to the target are removed from the data.It then shows phase-height versus time and phase-height rate estimates.A model is introduced to convert phase-height rate to AGB-rate, and AGB-rates are produced for the 78 plots.

Phase-Height Rate from Raw Interferometric SAR Data
The area of Tapajós analyzed (See Figure 1) was observed 32 times in 3.2 years (Table A1).The amplitude and phase were supplied by DLR for each spacecraft for small areas called "looks".A look extends 2.5 m in the azimuth direction (along the spacecraft flight line) and 1.4 m in the range direction (perpendicular to the flight line).Data presented in this format are called "single-look-complex".The data were referenced to radar coordinates of "azimuth" and "range", the number of looks along each direction.The amplitude-phase data, actually given as complex numbers with real and imaginary parts, were cross correlated by complex multiplication of the signals from each spacecraft e.g., [7] and averaged (as in ( 1)) over 50 × 50 m, 0.25-ha.In this process, an interferometric phase trend pertaining to the flat Earth was removed [28], leaving the phase of each 0.25 ha carrying only information about the altitude of the surface or vegetation.An amplitude and phase-height (phase/α z ) are available, for each of 32 epochs, for each of the 78 plots for which we have field-measured AGB.The amplitude and phase-height together have been used to estimate single-epoch AGB [10], but to estimate structural change in this work we use the phase-height only.As mentioned in Sections 2.2 and 2.3.3, each epoch of measurement corresponded to a different baseline, with a resulting different h amb , though the h amb 's are mostly in the 70-100 m range (Table A1).A reference epoch, called t 0 , of 22 September 2011 was adopted.Phase-heights from the reference epoch were subtracted from all other phase-heights and those differences, for each of 78 plots, were fit to time to extract a slope (rate).

Removal of Systematic Trends in Phase-Height Data
Before fitting phase-heights to time, a systematic, planar offset in phase-height for each baseline was removed.This offset depends on the range and azimuth, and has nothing to do with topography or vegetation structure, and therefore must be removed before measuring vegetation-related changes in phase-height.The offset is most likely due to unmodeled spacecraft and/or plot coordinates, as described below.This offset can be fit by a plane in the range and azimuth and then subtracted from the phase-heights [16,29].The planar model of the phase-height of plots i at time t j (h φ,i,j ) in terms of the phase-height of plots i at t 0 (h φ,i,0 ) is where range i , azimuth i are the radar coordinates of plots i, and a j , b j , and c j define the offset, range slope, and azimuth slope of the plane formed by the difference of phase-heights at t j and t 0 .They are estimated from the phase-heights of nearly all of the plots.The optimal solutions for 70 equations like (8) determine the plane parameters for each time t j .The jump stands, with abrupt discontinuities, are not used to estimate the parameters in (8), as they will cause rate biases, as explained in the next paragraph.
The overall offset a j results from one arbitrary phase offset per baseline, typical of interferometric data.The b j and c j slopes can result from unmodeled drifts in the spacecraft positions and/or imperfections in the flat-Earth phase removal mentioned above.In (8), δ i,j is the departure from the average plane of the phase-height for plots i at time t j .This residual is determined by inserting the a j , b j , and c j parameters estimated from 70 plots into an equation like (8) for each plot i.The phase-heights h φ,i,j and h φ,i,0 are known from the data, leaving the residual δ i,j the only unknown in (8).This residual for each stand, with one further correction, will be plotted versus time to determine phase-height rate.
The prescription of (8) yields phase-heights with the right spatial relative differences for a scene at a single epoch (a single j).For a series of times t j , it produces δ i,j 's with the correct absolute rate for plots i only if the rates averaged over all plots used to estimate a j were zero.For nonzero average rates over the 3.2 years of observation, the prescription of (8) yields relative rates of the δ i,j 's only.This is Forests 2017, 8, 277 11 of 28 because the time sequence of a j parameters will absorb the overall rate of growth of the plots used to estimate δ i,j , leaving the rate of the δ i,j 's artificially reduced (for positive overall growth) or augmented (for negative overall growth).In order to produce the absolute rates of the next section, two buildings were located in the InSAR data and their "growth rates" estimated as in the first line of (9) below from a time sequence of δ i,j 's from (8).The phase-height rate of the buildings should be zero, but the average of their phase-height rates was −0.45 m/yr.This negative phase-height change probably results from an overall growth of the vegetated plots used to do the estimation in (8).A constant 0.45 m/yr was added to the slopes of all time series of δ i,j by adding terms of the form 0.45 × (t j − t 0 ) to each stand's time series.The buildings' individual slopes were −0.4 and −0.5 m/yr, so 0.05 m/yr can be regarded as a possible systematic error in the rates of Figures 5 and 7.The means of arriving at this rate correction factor, by observing the time history of parts of the InSAR data known to be stationary (or some other known rate value), is potentially a way to assess overall change in the region spanned by the InSAR data, which appears to be positive (increasing vegetation height and consequent AGB; see Tables 1 and 2) for the sample of Tapajós in this study.This overall rate could in principle also be estimated by tagging trees and remeasuring heights for a few years, as described in Section 4.3.

Phase-Height Time Series and Rates
Figure 5 shows phase-height versus time for 9 examples drawn from the 78 field plots for which phase-height rate was estimated.Estimates of slope (phase-height rate) and slope error, along with the RMS scatter and AGB are shown in the text for each plot.The red lines show the best fit model, whether linear or the "jump" (i.e., logistic) model.Note that for some linearly fit plots, such as the one in the upper right or right middle of Figure 5, less than 3.2 years-1-2 years-of data would probably yield an acceptable rate error, of, say 0.5 m/yr.For most plots, a simple linear trend is the better fit, but for the top two middle plots in Figure 5, there was a clear jump.The two models, δ m,i,j used for plots i at time point j were δ m,i,j = int i + m i * t j linear and logistic function (9) where int i is the linear-fit intercept, m i is the phase-height rate parameter.For the logistic function, d i is a bias parameter, e i is the phase-height rate both before and after the jump in m/yr, f i is the "size" of the jump in meters, g i is related to the width or abruptness of the jump in yr −1 , and h i is a parameter specifying the epoch of the jump in years.The linear model is the default.The logistic model is chosen in Figure 5 if the size of the jump (the f i parameter) is greater than 4 m, and if the rms scatter of the logistic fit was at least 33% lower than that from the linear fit.A total of 8 disturbed plots were identified using these criteria.The fits adjusted parameters in (9) to minimize the reduced χ 2 using the errors bars in Figure 5, with reduced χ 2 for N observation epochs given by where σ i,j are the observational errors for stand i at time point j, and m is the number of parameters in the fit (2 for linear fits, 5 for logistic fits).Those observational errors were calculated by simulating the interferometric response of a forest of 100 randomly distributed particles and producing coherence and corresponding phase statistics.A lookup table was created giving phase-height scatter for each value of coherence, and phase-height error was thereby determined from the coherence data.For each stand, a single additional, unmodeled error, of order 0.5 to 1 m, was added in quadrature with the coherence-based error to each time point to bring the reduced χ 2 close to 1.The errors on the slope parameters shown for each plots are based on the observational phase-height errors.For the plots described by a linear fit, rate errors (after the ± symbols) were calculated using standard least-squares methods [30].For the plots described by the logistic function, errors in rates were calculated by Monte Carlo techniques.Errors in the h i parameters were about 1 month.Being able to identify the epoch of disturbance with this accuracy is a considerable gain for REDD+ monitoring over optical techniques which are plagued by nearly constant cloud cover in the tropics.Landsat, for example, had two clear observations in September 2013 and June 2014, compared to 11 over the same interval with TanDEM-X.Clearing indications were found in Landsat for all 8 plots with jumps.Some of the plots in Figure 5, for example the lower left, seem to show systematic effects which have few-month or annual time scales.Future analyses should include taking the Fourier transform of the traces in Figure 5 to see if annual or seasonal trends emerge.The rates for the complete set of 78 plots, along with their errors, latitudes and longitudes are shown in Table A2 in the Appendix A.
The rate of lidar RH90 was found to correlate with phase-height rate better than any other RH metric.The analysis which led to Figure 5 and Table A2 was redone including only those epochs between August of 2012 to August 2013 to match the lidar data.The TanDEM-X phase-height rate over this one year period is plotted versus the lidar RH90 rate in Figure 6  Figure 7 shows the estimated phase-height rates, a sample of which are portrayed in Figure 5, versus field AGB.The plots with AGB <250 Mg-ha −1 show a higher phase-height rate than the plots with AGB > 250 Mg-ha −1 .There seem to be significant peaks and dips for phase-height rates for AGB less than 100 Mg-ha −1 .A substantial number of points with AGB lower than 250 Mg-ha −1 are above the zero line by more than an error bar.Table 1 shows the average phase-height rate for secondary, primary, and all plots.The standard deviation and the standard deviation of the mean are shown in the 3 rd and 4 th columns.The standard deviation of the mean is the error in the average phase-height, given by the standard deviation divided by the square root of the number of plots in each category.Table 1 shows that secondary forest plots have a factor of ≈8 faster phase-height rate of change than primary plots.Some of the plots in Figure 5, such as the upper left panel, show statistically significant rates of change of phase-height, relative to the error in the rate parameter.Others, like the lower left panel, show insignificant rate.In Figure 7, many points are several error bars away from zero, and are therefore significant.Comparing columns 2 and 4 in Table 1, the average phase-height rate for secondary plots and that for all plots are significantly different from zero (column 2 is much greater than column 4).The average primary rates are consistent with zero.Note that systematic errors discussed in subsection 4.3 have not been addressed in the above discussion of significance.

From Phase-Height Rate to AGB Rate
This subsection derives the conversion from phase-height rate to AGB-rate.Given the closeness of phase-height to radar MCH discussed in Section 2.3, and given the many studies that have expressed AGB as a function of LiDAR MCH [20,[25][26][27], we start with AGB as an arbitrary function f of h φ and take the time derivative for A ĠB: Forests 2017, 8, 277 15 of 28 The conversion factor from phase-height rate to AGB-rate is the derivative of AGB with respect to phase-height.However, note that (10) refers to one plot.The derivative of AGB with respect to phase-height means the derivative with respect to changes in phase-height within a plot, not, for example, to changes in phase-height due to lateral displacement.As discussed in Section 3.4.1,this conversion factor can be taken as constant for each plot, over the ≈3-year period over which phase-height rates are estimated, for undisturbed plots.For jump plots, this approximation is strictly not true, in that AGB changes substantially for a jump and the conversion factor on the right side of jumps in Figure 5 should be different than that on the left.However, because the jumps happen in ≈2014, the linear rates are more established by the "before clearing" vegetation dynamics.In a more general formulation, the before-and after-jump AGB-rates could be calculated separately.In (10), the conversion factor is anticipated to be a function of AGB, as will be shown in Section 3.4.1.The conversion factor is derived in Section 3.4.1.The conversion factors are applied to the phase-height rate results to derive AGB-rates in Section 3.4.2.Section 3.4.3 is about the relative accuracy of single-epoch versus differential AGB estimates.

Model for Deriving Phase-Height-Rate to AGB Rate Conversion Factors
The conversion factor for transforming from phase-height rate to AGB-rate for a given plot is based on a model illustrated in Figure 8. Figure 8 shows an older (higher AGB) plot on the left and a lower-AGB plot on the right.They are both the same lateral size, 50 m × 50 m.The older one on the left has more growth than the one on the right, where the number of tree icons is schematically related to vegetation density at the top of the plot, z .The top height, z , is the height beyond which the mass density is zero.It is expected that the increment in AGB due to a change in top height z by ∆z will be larger for the left-hand plot than for the right-hand.This would certainly be true if, by increasing ∆z , the plot replicated another layer with the vegetation (and consequent mass) density characteristics at z .The increase in AGB would then be ρ(z )∆z , with ρ(z ) the plot's mass density at z .However, it is also possible that the increase in mass density in going from z to z + ∆z occurs lower in the canopy, near the base of the trunks, for example.In general, therefore, ρ(z ) is to be interpreted as the quantity by which ∆z must be multiplied to give the mass/area increase of the plot when the top altitude z → z + ∆z .
Δ"′ Δ"′ z' z' In either case, the AGB needed to evaluate the derivative in ( 10) is an integral over top heights z : where ρ(z ) is now restricted in the second line of (11) to be a power law for each stand.The rationale for plot-level mass density being proportional to a height metric to a power is loosely based on the notion that individual tree mass goes as a power of height [23].However, that a collection of power-law trees sum to a power-law plot has not been demonstrated; it is taken as an empirical model.The exponent β in (11) in the literature has taken on values of 1-3 [12,14,16], exponentiating total height, Lorey's height, or phase-height.It is conceivable that β varies from plot to plot, or has an AGB dependence.In what follows, however, β will be taken to be common to all plots.The amplitude is η, which will be allowed to vary between plots.In (11), h v is the total physical height.We will express the conversion factor in terms of phase-height below.
In order to evaluate the conversion factor derivative in (10), a small change in a single plot's AGB, ∆AGB is considered in response to a small change in total height, ∆h v : The motive for writing the conversion factor in curly brackets as in the last line of ( 12) is to express A ĠB in terms of accessible quantities from TanDEM-X.Because it is expected (and suggested in Section 4.2) that model-based total physical height estimates perform poorly for AGB-rate estimation from this data set, we use phase-height, a more direct product of TanDEM-X.This study is therefore pitched at potential missions where phase-height is readily available, while multiple baselines and polarizations needed for total-height estimation may not be available.The term in square brackets will be applied as an over-all correction factor.For no attenuation or infinite attenuation, the term in square brackets is unity.A random-volume model calculation shows that the term in square brackets depends on extinction and forest height, neither of which are readily available for all of our plots via remote sensing.Averaging over extinction coefficients from 0.1 to 0.5 dB-m −1 and heights from 0 to 40 m, it was found that the term in square brackets is on average 0.85 with a standard deviation of 0.05.In the absence of detailed extinction or h v information, this 0.85 factor was applied to all stands and ±0.05 will be considered a source of error.
From (12), in order to arrive at the conversion factor, the quantity AGB/h φ must be estimated.AGBs from field measurements were combined with two sources of h φ .In previous work, 30 plots from Tapajós (our first 26 + 4 more outside of the rectangle of Figure 1) were analyzed for h φ .It was obtained visually by looking at the lower extrema of look phase-heights to calibrate the ground phase-height [10].The ratio of the field-measured AGB to phase-height for those plots is plotted in Figure 9 in red.Another potential source of phase-height can be gotten from the jump plots in this analysis.For each of the 8 plots with jumps, the phase-height change can be considered the total phase-height before the clearing event; this is true if all of the AGB was taken away.AGB/h φ ratios from those 8 plots are shown in black in Figure 9.  12) AGB/phase-height.This term is proportional to the derivative of AGB with respect to phase-height, if a polynomial model is used as in the second line of (11).The red points are from InSAR phase-heights calibrated with visual ground finding in another study [10].The black points are the " f i " parameters (see ( 9)) from logistic-fit jump events.The dashed red line is the fit in ( 13) based on the red calibrated points.Because the black points seem to lie consistently above the red points, it was assumed that the total AGB was not removed for some of them.They were not used in developing the conversion factor in (12).
The black points being consistently higher than the red points calibrated for ground location suggests the possibility that all the AGB was not removed in some of the jump events.It was determined that about 25% of the AGB would have to have been left behind to get the black points to lie on the red ones in Figure 9.The red, calibrated points were therefore chosen for the dashed-line fit, which describes the AGB/h φ term in the conversion factor as a function of AGB: The RMS AGB/phase-height about the line in (12) was 0.3 Mg-ha −1 -m −1 and is considered a source of error in Section 4.3.This small RMS suggests that the correction factor (β multiplied by (13) times the square bracket term in (12) of 0.85) can be considered just a function of AGB and need not be considered a function of phase-height.This is probably because, for any AGB, the range of possible h φ s is restricted.Thus, specifying AGB automatically specifies the allowed h φ values.
AGB could be provided by, for example, the Global Ecosystem Dynamics Investigation (GEDI), but acquiring spatially and temporally dense phase-height from any other sensor but a radar interferometer would be much more difficult operationally.That the conversion factor can be modeled as a function of AGB only could have beneficial operational consequences.However, measurements of AGBs and h φ s would have to establish a curve like (13) for each region.The universality of (13) is not known; a plot like Figure 9 should be generated from other tropical areas.The conversion factor at Tapajós ranged from 1 to 9.5 Mg-ha −1 -yr −1 /(m-yr −1 ) for the jump plots, and up to ≈17 for the calibrated plots.
For each plot, the conversion factor can be considered constant over the 3.2-yr observation period, if the AGB in the conversion factor is taken to be that at the middle of the observation period (3.2 years in our case).If, for example, a 10 Mg-ha −1 plot went from 0 Mg-ha −1 to 20 Mg-ha −1 over 3.2 years, the asymmetry of ( 13) about the midpoint of 10 Mg-ha −1 would cause 1% error in A ĠB.For higher AGBs, the asymmetry effect would be even smaller.

Applying Conversion Factors to Estimate AGB Rate
Figure 10 shows AGB-rate, calculated from the points in Figure 7, with conversion factors based on AGB as in ( 12) and ( 13), with β = 1.We take β = 1 in part because it has been used in the literature [14,16], but mainly because the results are then easily scalable to any value of β.The results using β between 1 and 2, as will be seen below, are within a reasonable range of other external measures of tropical-forest rate.From ( 12), the conversion factor for an arbitrary β is just the product of that β and the β = 1 factor.Using β = 1, the average rate is 1.7 Mg-ha −1 -yr −1 with a standard deviation of 4.0 Mg-ha −1 -yr −1 .The distribution of rates is shown in Figure 11.The AGB rates in Figure 10 for the complete set of 78 plots, along with their errors are shown in Table A3 in the Appendix A. The AGB rates estimated from the phase-height rates multiplied by the conversion factor in (12), plotted against AGB.The green curve is the site-specific, stand-level growth model of Neeff and Santos [31] based on field data at Tapajoś.

Applying Conversion Factors to Estimate AGB Rate
Figure 10 shows AGB-rate, calculated from the points in Figure 7, with conversion factors based on AGB as in ( 12) and ( 13), with β = 1.We take β = 1 in part because it has been used in the literature [14,16], but mainly because the results are then easily scalable to any value of β.The results using β between 1 and 2, as will be seen below, are within a reasonable range of other external measures of tropical-forest rate.From ( 12), the conversion factor for an arbitrary β is just the product of that β and the β = 1 factor.Using β = 1, the average rate is 1.7 Mg-ha −1 -yr −1 with a standard deviation of 4.0 Mg-ha −1 -yr −1 .The distribution of rates is shown in Figure 11.The AGB rates in Figure 10 for the complete set of 78 plots, along with their errors are shown in Table A3 in the appendix.The AGB-rates for AGB < 250 Mg-ha −1 are on average positive, and become close to zero for AGB > 250 Mg-ha −1 .Table 2 below summarizes secondary-, primary-, and all-forest average AGB-rates.The AGB-rates for AGB < 250 Mg-ha −1 are on average positive, and become close to zero for AGB > 250 Mg-ha −1 .Table 2 below summarizes secondary-, primary-, and all-forest average AGB-rates.Using the same considerations as after Table 1, the secondary-forest and all-plot average AGB rates are significantly different from zero; the primary-forest AGB rate is not.
It should also be noted that if a curve were fit through the jump plots, the black points in Figure 9 instead of the red, average conversion factors and AGB-rates could change by ≈ 30%.Biases in ground finding in producing the red points could contribute to the red points being lower than the black.The selection of red or black points to produce Figures 10 and 11 therefore is a potential source of systematic error at the level of 30%.
The green curve in Figure 10 is the stand-level growth model of Neeff and Santos [31], which was based on field measurements at Tapajós and allometrics.It was derived from chronosequences.The curve demonstrates that the AGB-rates estimated here are of the order of the modeled rates, with the model giving higher growth rates for the secondary plots, possibly arguing for β > 1.The reference does not supply rates for AGB > 200 Mg-ha −1 .Note that Figure 10 is also consistent with Feldpausch [32] where upper limits of about 10 Mg-ha −1 -yr −1 are quoted for tropical moist forests.These four external measures could be seen as constraining the value of the exponent β in the second line of (11) to between 1 and 2. Both the TanDEM-X AGB-rates and those of the green-line model show a rising rate from 0 to about 80 Mg-ha −1 -yr −1 .The model takes a downturn at ≈80 Mg-ha −1 , the TanDEM-X AGB-rates show some indication of a downturn as well, but then show a higher rate contribution starting at around 200 Mg-ha −1 , and eventually decreasing for more massive plots.The error bars are the slope errors multiplied by the correction factors of (12).Recall that the slope error bars include an added unmodeled error to make the reduced χ 2 close to 1.The AGB-rate error bars therefore reflect the goodness of fit (RMS scatter) of the linear or logistic model.There appear to be significant trends, given the error bars.For example, as noted above, around 200 Mg-ha −1 there appears to be a collection of points significantly positive (their error bars do not clip zero).The significant departures of points from any group could be 1) due to using different ambiguity heights for different epochs, as mentioned in Section 2.3.3.While the error introduced by using different ambiguity heights seems small, based on considerations in Section 2.3.3, it is possible that pathological asymmetries in radar power profiles, p(z ), could couple with multiple h amb s to cause experimental errors.In addition to experimental errors, some of the apparently significant points in Figure 10 could in fact 2) be indicative of real structual dynamics.They could also result from 3) seasonal dielectric change, which could change p(z ) and thereby phase-height.In the absence of plot-by-plot field rate measurements, we cannot discriminate against these three possibilities.

The Relative Accuracy of Single-Epoch and Differential (Rate) AGB Estimates
Each point in Figure 10 is derived from a phase-height versus time graph, as in Figure 5, converted to AGB versus time by multiplying by the conversion factor in (12).The RMS scatter of the AGBs in a fit is the RMS of the phase-height fit multiplied by the conversion factor.The average AGB RMS in those 78 fits is 12 and 24 Mg-ha −1 , for β = 1 and β = 2, respectively.In contrast, using phase-height and coherence to estimate single-epoch AGB, the RMS of 30 plots' AGB about field-estimated AGB is 57 Mg-ha −1 [10] with data from the same Tapajós set.If the linear (or logistic) trend here and the field-measured AGB in the other work are regarded as truth, the scatter about the truth is roughly 2-4 times greater for the single-epoch measurement than for the rate (differential) measurement.This could be due to the fact that the location of the ground had to be estimated for the single-epoch measurement for each stand.For the differential measurement, the difference between each phase-height and that of the reference date, ground-finding errors common-mode cancel.In addition to effects of imperfect ground finding, single-epoch studies are subject to plot-to-plot differences in complex structure.The radar response to those differences-whether in coherence, phase-height, or power-can trigger erroneous biases in AGB estimated from the radar signals at a single epoch.In contrast, the differential measurement is not sensitive to plot-to-plot structural biases, as long as those biases persist over time.

Discussion
The methods and results of this paper suggest that measurements of X-band phase-height are potentially an indicator of tropical-forest structure and AGB change.The efficacy of using X-band in tropical forests with interferometry will be discussed below, as well as the use of phase-height from InSAR.Future improvements to the study described here follow.

X-Band Interferometry for Tropical-Forest Structure and Biomass
It has long been thought that X-band signals would not penetrate the tropical forest sufficiently to sense AGB much above 30 Mg-ha −1 [33].Because interferometry involves phase, it has been suggested that interferometry, even at higher frequencies, has more of a "height perception" [34].A simple way to envision the gains by using phase is to imagine that the forest is a nearly impenetrable layer.If radar only sees the top few cm, say, power measurements will be the same if that impenetrable layer is displaced vertically by a few meters.However, if the nearly impenetrable layer is raised by 1 m, the InSAR phase-height will change by 1 m.InSAR is sensitive to the spatial origins of the return radiation, and this is the big difference between interferometry and traditional power-based radar.The origins of return radiation-whether from InSAR or LiDAR [35]-bear on AGB estimation.In fact, the X-band signal does appear to penetrate, yielding a vertical diversity of returns.Although not used in this paper, X-band coherence should decrease as plot AGB goes up, if there is penetration.Coherence was shown to decrease with increasing AGB for a subset of the plots in this study [10].
Signals attenuate more at higher frequencies principally because the absorption by water increases with frequency [36].However, as frequency goes up, the ability of the signal to penetrate holes in the medium increases.Thus modeling suggests that an X-band signal appears to be able to penetrate 50 cm holes e.g., [37], even though it may be severely attenuated by vegetation in its path.The inspection of TanDEM-X data suggests that the tropical-forest signal at X-band is consistent with hard, attenuating targets and significant gaps, allowing more penetration than previously suspected [38,39].

Using Phase-Height to Measure Tropical Forest Dynamics
The time evolution of phase-height was chosen for AGB-rate estimation.It was chosen because phase-height is close to the radar-power-averaged mean canopy height, which, in turn, is plausibly connected to profiles of the density of scatterers, as detailed in Section 2.3.Phase-height was chosen in part because the median 1.3-m RMS scatter about a linear or logistic model is smaller than errors attained using other height metrics with InSAR alone, as detailed below.
Phase-height is one of two observations from interfermetric SAR.Coherence, the normalized amplitude of the InSAR cross correlation in (1), is the other.Coherence has been used along with phase-height in simple models-e.g., Random Volume Over Ground (RVOG)-to estimate forest height, extinction, and underlying topography [7,40,41].Forest height RMS's of RVOG height about LiDAR had values of 2 m in boreal forests, 2-4 m in tropical forests when used in conjunction with LiDAR, and >5 m for InSAR alone in tropical forests [42].In [17], 3-4 m RMS's for forest height were demonstrated with (3.5 m) and without (4.3m) LiDAR in temperate forests.TanDEM-X phase-height, in addition to its correspondence with MCH and density profiles, produces lower RMS's than most of the height estimates arrived at by modeling.This study investigated a monitoring scenario with the temporal frequency of Figure 5, and therefore focused on "InSAR-only" phase-height rate monitoring, which motivated the focus on the high-performing phase-height for dense temporal measurements.This work initiates a thorough understanding of phase-height for change measurements, and begins to ask whether that might lead to approaches for monitoring dynamics that are complementary to model approaches such as RVOG.The RVOG performance might be improved if there were more frequent HV (horizontal send, vertical receive) polarization data, and simultaneous, baseline acquisitions.The RVOG attempts to date with TanDEM-X have used only the dominant HH and VV signals at one baseline per epoch, and this could be part of the reason for the poorer performance.It could also be that the RVOG model best estimates structure for more penetrating frequencies, such as L-band, for which ground returns will contribute more polarimetric diversity, and thereby lower the errors on height estimates.

The Performance of Phase-Height and AGB Rate and Future Enhancements
All estimates of the performance of phase-height and AGB rates for measuring slow change are based on a 3.2 year observing period, with 32 observations, as in Figure 5.An indication of the phase-height rate performance is the array of formal error bars in Figure 7. Recall that these errors result from phase-height errors, which had an added component to account for unmodeled scatter.The average formal error is 0.27 m-yr −1 with a standard deviation of 0.15 m-yr −1 .The scatter of phase-height rate about the y = x line in Figure 6 provides another approach to assessing phase-height rate error.The 1.2-m-yr −1 scatter of InSAR phase-height rates about LiDAR RH90 rates can be interpreted as the 1-yr error in phase-height rate (assuming LiDAR errors are much smaller than those of InSAR).Assuming a time interval −1.5 dependence of the rate error estimate, an error estimate of 0.21 m-yr −1 is derived for a 3.2-yr interval, which is of the order of the formal error bar based on the phase-height fit alone.Similarly, the performance of AGB-rates from the array of error bars in Figure 10 is indicated by the average error of 1.9 Mg-ha −1 -yr −1 with a standard deviation of 2.2 Mg-ha −1 -yr −1 .From an average phase-height rate error of 0.27 m-yr −1 , with an average conversion factor of 7, an AGB-rate error estimate of about 1.9 Mg-ha −1 -yr −1 results, equal to the rate error above.A source of systematic error in the phase-height rate is in the correction derived from the two buildings.The correction of 0.45 m-yr −1 was the average of 0.5 and 0.4 m-yr −1 of the two buildings.A phase-height systematic error of about 0.05 m-yr −1 can be assumed in the phase-height rate, which, multiplied by the average conversion factor of 7, gives a possible AGB-rate systematic error of ≈0.4 Mg-ha −1 -yr −1 .In the future, many buildings or ground control points could be sought to drive this error lower.
The AGB-rate error of 1.9 Mg-ha −1 -yr −1 is based on the scatter about a linear or logistic phase-height fit and does not account for important systematic errors in the conversion from phase-height rate, including model assumptions in ( 10)- (12).The largest systematic effect in AGB-rate is the choice of β in (11).The green curve in Figure 10 argues for a higher β to bring the low-AGB red points closer.The study which generated the green line did not use the same plots as in this study, nor was it done at the same general epoch.Other studies mentioned in Section 3.4.2suggest limiting β.A range of 1-2 seems plausible.The black curve of Figure 9, based on cleared areas in this study, suggests a 30% increase in the conversion factor, and hence all AGB-rates.In general, there are many alternatives to the treatment of the conversion factor in (10).The density function, for example, in (11) could be some other function of z and not a power law at all.Differential fieldwork, on so-called "permanent plots" with tagged trees, could suggest other functional forms for ρ(z ), including values for β.The scatter of the red points about the dashed line of Figure 9 is 0.3 Mg-ha −1 -m −1 , which, when multiplied by the average phase-height rate (from Table 1) of 0.5 m-yr −1 yields a 0.15 Mg-ha −1 -yr −1 systematic effect.This scatter is in part due to the error in measured AGB of 25% [24].In looking up a plot's AGB/h φ , a 25% error will be made in the AGB argument of (13).This translates to a 15-25% error in AGB/h φ .
Another source of systematic error is the treatment of the square-bracket term in (12).It was assigned an average value of 0.85, with a standard deviation of 0.05, based on modeling attenuation and a random volume.Using this value of 0.85 therefore constitutes a 6% systematic error.The permanent-plot fieldwork could help to model this term as a function of AGB, and apply the functional form of the term without having to resort to just an average value.
Most systematic errors would take the form of overall scale factors.Given that all the systematic effects (except for assignment of β) cause <1 Mg-ha −1 -yr −1 errors, the AGB-rate error from the error bars or LiDAR comparison of ≈2 Mg-ha −1 -yr −1 quoted here is plausible, but probably an underestimate by 1-2 Mg-ha −1 -yr −1 .
The best indicator of the error is a comparison to field estimates of change, as mentioned repeatedly.In fieldwork done in 2010 and 2013, there were a few plots repeated, but not tagged, so AGB errors were about 25%, as noted above, and therefore too large to be an indicator of truth AGB-rate.We never intended to re-measure the 2010 plots and that is why they were not marked.The comparison to external reports mentioned just below Table 2 shows that all approaches, in Tapajós, Costa Rica, and Panama, get the same order of secondary-forest and primary-forest change [19,20].
Estimates of the performance of phase-height for measuring jumps, or abrupt disturbance, entail errors in the epoch of the event and the magnitude of the phase-height change.From Figure 5 and Monte Carlo error analysis, the epoch and magnitude of disturbance can be determined with about 1 month and 2 m accuracy, respectively.All of the 8 jump events, 2 of which are shown in Figure 5, were seen in Landsat images, with temporal resolution of about 1 year.The ability to establish a 1-month window for clearing events seems a strength of phase-height measurements, and will be of considerable use to REDD+ mitigation activities.
On the data acquisition side, an improvement could be to use the same h amb for each epoch.The model calculations leading to Figure 4 suggest that even with extremely asymmetric radar power profiles, only a 0.6 m difference arises between MCH and InSAR phase-height.This 0.6 m could be contributing to the RMS phase-height scatters in 5, but it seems unlikely as the MCH-h φ comparison was a worst-case result for the volume calculation of (7).It is also possible that pathological asymmetries could result from adding ground contributions to the model in (7).In one variant of the analysis, extremely large height ambiguities of 250 m were removed without appreciable difference in phase-height rate.It does appear that phase-heights from similar h amb exhibited some clumping in plots like those of Figure 5. If, in the future, all h amb s were the same, this worst-case RMS contribution would drop to 0.1 m for the phase-height scatter.
Although the revisit intervals and accuracy of current TanDEM-X measurements is sufficient to detect the epochs of observed jumps to within about a month, the error estimates above for measuring growth or gradual degradation of 2-4 Mg-ha −1 -yr −1 or more apply to a 3.2 year period.The average ratio of the magnitude of estimated AGB slope to slope error bar is 2.9.Assuming a time-span −1.5 error dependence again, almost a factor of two accuracy would be lost in going to a 2-yr delivery period.This factor of two could probably be recouped by using 1 ha sample instead of 0.25 ha of this study.

Conclusions
This study suggests that InSAR phase-height rate from TanDEM-X, which is close to radar MCH rate, can be determined by time series of phase-heights to about 0.25 m-yr −1 .In Tapajós, secondary plots averaged a phase-height rate of 0.8 m-yr −1 , while primary plots averaged 0.1 m-yr −1 .This study further proposes a model in which AGB-rate is proportional to phase-height rate, with the conversion factor between the two rates dependent on AGB.AGB-rate errors of 2-4 Mg-ha −1 -yr −1 for 3.2-year time spans and 0.25-ha plots were based on error modeling.The average secondary-forest AGB-rate was 2.1 Mg-ha −1 -yr −1 and primary-stand rate of 1.1 Mg-ha −1 -yr −1 .The AGB-rates are subject to as much as a factor of 2 systematic effects due to power-law model assumptions and other systematic effects as in section 4.3.Other systematic effects are less than 1 Mg-ha −1 -yr −1 .The epochs and sizes of disturbance events were determined with 1-month and 2-m accuracy, respectively.Permanent plot fieldwork will facilitate analysis approaches with smaller systematic effects.TanDEM-X is scheduled to provide more forest data in 2017.If this continues to 2018 and beyond, it affords an opportunity to do accurate permanent-plot fieldwork to support new TanDEM-X acquisitions.
The scenario suggested in this paper for converting phase-height rate to AGB-rate (10) requires a single-epoch measurement of AGB for each plot as input to the conversion factor in (12).More elaborate models than (10)- (12) will probably also require a single-epoch measurement.The InSAR coherence and phase-height from TanDEM-X or a similar mission could be used to generate a 30% AGB as in [10].However, in that work, a by-hand method was used to find the ground to calibrate the phase-height (to realize the red points in Figure 9).An automated method, sufficiently accurate, would have to be developed to facilitate the AGB estimates.Another possibility is to use single-epoch LiDAR measurements of AGB from a GEDI-type mission, along with phase-height from an X-band InSAR mission, to arrive at conversion factor components as in (13).GEDI measurements in the tropics will probably be infrequent due to cloud coverage, but for undisturbed plots, an accurate measurement of AGB within 1 year will be sufficient (AGB will not have changed, on average, by more than 2-3 Mg-ha) to find the conversion factor.Combining LiDAR's accurate but infrequent AGB measurement, with InSAR's frequent phase-height, could be an efficient scenario.Yet another possibility for accurate single-epoch measurements would be an L-band, fully polarimetric 2-spacecraft InSAR mission, which could undoubtedly find the ground better than an X-band mission.
The approaches presented here, and alternatives to them which might use model-based height metrics from InSAR, should be tried at various radio frequencies.There is no firm argument for a correspondence between the "hard scattering" at X-band discussed in Section 4.1 and the performance of phase-height rate and AGB-rate.The estimated heights from InSAR based on RVOG perform poorly at X-band with low polarimetric and baseline diversity, as noted in Section 4.2.However, heights from RVOG or other models will almost certainly produce more accurate single-epoch stock measurements at lower radio frequencies such as L-band.Whether the meter-level scatters in Figure 5 at X-band would improve at L-band should be studied.More elaborate models of the dependence of A ĠB on ḣφ than (10)

Figure 1 .
Figure 1.The Tapajós National Forest covered by TanDEM-X (large white rectangle) in this study.Yellow dots are 78 field plots, and the smaller red rectangles are areas of LiDAR coverage.

2 Figure 2 .
Figure 2. Schematic representation of the InSAR scattering geometry.The fields scattering from a single scattering element, the top of the tree, travel along the distances r 1 and r 2 to the ends of the interferometer.This difference is proportional to interferometric phase and to the height of the single scattering element.The interferometric phase resulting from many scattering elements (the volume of the forest) is used to characterize vertical structure in this paper.

2. 3 . 3 .
Quantitative Relationship: The Magnitude of Differences Between Phase-Height and Mean Canopy Height Figure 4  was generated by starting with σ right = 2 m, which is the symmetric profile of p(z ) versus altitude in meters, in the lefthand inset.The right-hand standard deviation, σ right , was then stepped from 2 m to 12 m, which resulted in the maximally asymmetric profile in the right-hand inset.Any further extension of σ right would lead to a more symmetric profile.As σ right was increased from 2 m, the MCH increased and was plotted on the x-axis.The y-axis is the phase-height-MCH difference obtained by inserting (7) into (3) and(2).The resulting h φ -MCH differences are shown for two values of h amb , 60 m and 250 m, which nearly bracket the range of h amb used in the h φ -versus-time results in the next section (see Table

Figure 4 .
Figure 4.The calculated difference between phase-height and mean canopy height (power-profile-averaged height), realized for an asymmetric-Gaussian hypothetical model of the volume radar power profile p(z ) as in(7).σ le f t is held at 2 m.When σ le f t = σ right , p(z ) is symmetric, the h φ -MCH difference is zero, and MCH = 12.5 = z 0 .The symmetric profile of p(z ) versus altitude in meters is in the left-hand inset.With σ right = 12, the right-hand inset shows the highly asymmetric power profile as a function of altitude in meters, and the h φ -MCH difference is −0.6 m.The h φ -MCH difference is plotted against MCH, which goes from 12.5 to 20 m when σ right is increased from 2 to 12 m.The difference is shown for two heights of ambiguity.

Figure 5 .
Figure 5. Phase-height versus time in years for 9 of 78 plots.All heights are differences relative to 22 September 2011.Red lines are linear fits of phase-height to time for plots without a jump."Slope"is the fit linear rate.The number after the ± sign is the formal error on the slope.Plots with mass <250 Mg-ha −1 tend to have more positive rates than those with AGB > 250 Mg-ha −1 .Two plots with jumps, due to agricultural clearing are in panels 2 and 5.The red line for those plots is the logistic function.For those plots, "slope" is a single rate estimated with the logistic function applied before and after the jump.The AGB (before clearing for jump plots) is also shown, and the RMS about the fit function is also shown.
, demonstrating a 0.5 correlation, with a probability of 1 part in 10 4 of realizing the Figure by accident.The RMS about the y = x line is 1.2 m-yr −1 .

Figure 6 .
Figure 6.TanDEM-X phase-height rates, evaluated over August 2012-August 2013 period versus RH90 rates estimated by lidar (red rectangles in Figure 1) over the same time period.The dashed line is y = x.

Figure 7 .
Figure 7.The phase-height rates estimated from (9) versus field AGB.Many of the points below 250 Mg-ha −1 are above the dashed zero line by more than an error bar, suggesting higher phase rates at lower values of AGB.

Figure 8 .
Figure 8. Schematic illustration of two stands with an increment of ∆z in top height.The number of tree icons represents vegetation and consequent mass density in the stand near the physical top height z .The same change in top height, ∆z , will produce less AGB change in the shorter plot on the right than in the taller one on the left.

Figure 9 .
Figure 9.The term in the conversion factor (12) AGB/phase-height.This term is proportional to the derivative of AGB with respect to phase-height, if a polynomial model is used as in the second line of(11).The red points are from InSAR phase-heights calibrated with visual ground finding in another study[10].The black points are the " f i " parameters (see (9)) from logistic-fit jump events.The dashed red line is the fit in (13) based on the red calibrated points.Because the black points seem to lie consistently above the red points, it was assumed that the total AGB was not removed for some of them.They were not used in developing the conversion factor in(12).

Figure 10 .
Figure10.The AGB rates estimated from the phase-height rates multiplied by the conversion factor in(12), plotted against AGB.The green curve is the site-specific, stand-level growth model of Neeff and Santos[31] based on field data at Tapajoś.

Figure 10 .
Figure10.The AGB rates estimated from the phase-height rates multiplied by the conversion factor in(12), plotted against AGB.The green curve is the site-specific, stand-level growth model of Neeff and Santos[31] based on field data at Tapajoś.

Figure 11 .
Figure 11.Histogram of the AGB rates of Figure 10 for the 78 plots in this study.

Figure 11 .
Figure 11.Histogram of the AGB rates of Figure 10 for the 78 plots in this study.

Table 1 .
Average, deviation, and standard deviation of the mean of phase-height rate for secondary forests, primary forests, and all plots.

Table 2 .
Average, standard deviation, and standard deviation of the mean of AGB-rate for secondary, primary, and all forest plots.