An Automatic Mosaicking Algorithm for the Generation of a Large-Scale Forest Height Map Using Spaceborne Repeat-Pass InSAR Correlation Magnitude

Yang Lei and Paul Siqueira *Department of Electrical and Computer Engineering, University of Massachusetts Amherst,151 Holdsworth Way, Amherst, MA 01003, USA; E-Mail: ylei@ecs.umass.edu* Author to whom correspondence should be addressed; E-Mail: siqueira@ecs.umass.edu;Tel.: +1-413-577-0623; Fax: +1-413-545-4652.Academic Editors: Nicolas Baghdadi and Prasad S. ThenkabailReceived: 5 February 2015 / Accepted: 13 April 2015 / Published: 5 May 2015Abstract: This paper describes an automatic mosaicking algorithm for creating large-scalemosaic maps of forest height. In contrast to existing mosaicking approaches through usingSAR backscatter power and/or InSAR phase, this paper utilizes the forest height estimatesthat are inverted from spaceborne repeat-pass cross-pol InSAR correlation magnitude.By using repeat-pass InSAR correlation measurements that are dominated by temporaldecorrelation, it has been shown that a simpliﬁed inversion approach can be utilized to createa height-sensitive measure over the whole interferometric scene, where two scene-wideﬁtting parameters are able to characterize the mean behavior of the random motion anddielectric changes of the volume scatterers within the scene. In order to combine thesesingle-scene results into a mosaic, a matrix formulation is used with nonlinear least squaresand observations in adjacent-scene overlap areas to create a self-consistent estimate offorest height over the larger region. This automated mosaicking method has the beneﬁtof suppressing the global ﬁtting error and, thus, mitigating the “wallpapering” problem inthe manual mosaicking process. The algorithm is validated over the U.S. state of Maine byusing InSAR correlation magnitude data from ALOS/PALSAR and comparing the invertedforest height with Laser Vegetation Imaging Sensor (LVIS) height and National Biomassand Carbon Dataset (NBCD) basal area weighted (BAW) height. This paper serves as acompanion work to previously demonstrated results, the combination of which is meant to


Introduction
Large-scale mosaics of biophysical parameters (such as biomass and forest height) are important for monitoring global carbon storage, as well as forest degradation.For this purpose, spaceborne missions have the vantage point from space and are efficient platforms for remote sensing data collection.The Shuttle Radar Topography Mission (SRTM), which carried a C-band synthetic aperture radar (SAR) with an interferometric (InSAR) baseline of 60 m, was on-board the Space Shuttle Endeavor during an 11-day mission in the year 2000.The data products from SRTM are processed into global mosaic maps of digital surface models (DSM) [1].When combined with field inventory and optical data, the SRTM InSAR phase data has been utilized to generate national mosaic maps of both biomass and forest height for the United States, which are provided by the National Biomass and Carbon Dataset (NBCD) 2000 [2].In addition to this data, JAXA's L-band JERS-1 (from 1992 to 1998; [3]) and ALOS/PALSAR (2006 to 2011; [4,5]) have collected a large amount of SAR image data over the past decade.Continental-scale mosaics of SAR backscatter intensity have been demonstrated and related to both spatial and temporal forest characteristics [6][7][8].Mosaics of polarimetric SAR backscattering coefficients can either be converted to forest biomass maps through a regression method [9,10] or classified into forest/non-forest maps [7,11], where the detailed forest structural characteristics can be further obtained through image segmentation techniques [12,13].
A simple mosaicking method for SAR/InSAR observations makes use of ephemeris information from the satellite and is a scene-by-scene approach that uses scene-overlap regions to estimate positional errors in the scene locations.When implemented using a limited number of scenes as the reference location, errors in the geometric locations are propagated away from the references as more scenes are added [14].The errors in the geometric offsets of individual scenes can be considerably and simultaneously reduced through the use of overlap areas between adjacent scenes when a matrix formulation of the geometric transformation is established [1,6].Terrain information from accurate DSM data (e.g., SRTM) is often used in this context to correct both geometric and radiometric errors in the data [7,8].
In contrast, this paper deals with a slightly different mosaicking problem, which comes from a simplified forest height inversion model [15].This model works through utilizing the spaceborne repeat-pass InSAR correlation magnitude (neither InSAR phase nor SAR backscatter power, as introduced above) measurements that are dominated by temporal decorrelation effects.Although spaceborne single-pass InSAR systems are capable of extracting forest height information [16,17] from volume scattering effects without encountering the dominant error of temporal decorrelation, results from such systems have not been demonstrated over large areas.
In order to make use of available repeat-pass InSAR measurements, where the temporal decorrelation (due to both dielectric change and random motion) is the dominant factor, a repeat-pass InSAR scattering model can be used to estimate forest height from this temporal decorrelation [15].In this model, the temporal change effects from both the dielectric change and the random motion among scatterers are incorporated into the random volume over ground (RVoG) model [16,17], which is adapted to repeat-pass InSAR correlation measurements.The dielectric change effect is included as a height-independent term (and, thus, independent of the volume scattering), while the random motion effect has a Gaussian height-dependent profile, and is found to be tightly coupled with the volume scattering.Therefore, as long as there is a vegetation-height dependence for temporal decorrelation, it is possible to utilize this temporal decorrelation to invert forest height, even without volume decorrelation, which is accomplished by assuming that: (1) the temporal change effects and forest backscatter profile/extinction coefficient follow some mean behavior across each interferogram; (2) there is minimal ground scattering contribution for HV-polarization; and (3) the interferometric vertical wavenumber is small.Under these assumptions, it is noticed that the temporal change effects dominate, and therefore, a simplified inversion approach is developed to link the observed HV-polarized InSAR correlation magnitude to forest height through the use of some known ground validation height data (e.g., LiDAR).The introduced model parameters of the inversion approach are: S scene (0 ≤ S scene ≤ 1), which characterizes the correlation component of dielectric change, and C scene (usually 0 ≤ C scene ≤ 20 for ALOS scenes), which primarily describes the level of random motion.The model parameters derived from this supervised regression at the ground validation site consisting of 44,000 hectares in central Maine are used as the basis for propagating the estimates of forest height to available interferometric pairs (through the overlaps between adjacent InSAR scenes) for the entire state (about nine million hectares), thus creating a state-mosaic map of forest height.
The mosaicking procedure used in [15] is outlined in Figure 1 for 37 ALOS/PALSAR scenes spanning over the entire state of Maine.Each scene is named by its ALOS orbit number followed by its frame number, with the color implying the stage of the above-mentioned manual mosaicking process, while each directed edge shows the direction of propagation through the overlap areas.There exists a small strip of Laser Vegetation Imaging Sensor (LVIS) [18] data over the Howland Forest in the central Maine area (marked in "red" and named "119_890").Data from the overlap between the InSAR height model and the LVIS-measured height are then manually propagated by sequentially going through the adjacent-scene overlap areas, which leads to the mosaic map of forest height for the entire state of Maine.
This manual mosaicking approach is prone, however, to what can be termed as the "wallpapering" problem (that is, by fixing one or two points of the wallpaper and gradually attaching the remaining part, larger deviation will occur as distance from the fixed points increases).Here, scenes that are farther away from the ground validation sites result in larger uncertainty in the determined estimates of the model parameters.In addition to this effect, the propagation path/sequence is non-unique, leading to a non-unique solution in determining the model parameters.The solution to this problem is to introduce an automatic mosaicking algorithm that estimates the desired model parameters simultaneously and arrives at a solution that is mathematically traceable and has a globally-minimized error.Figure 1.The outline of the mosaicking scheme.Each interferogram is represented by "orbit#_frame#".The mosaicking process starts from "119_890" (marked in "red") in the central Maine area, where Laser Vegetation Imaging Sensor (LVIS) LiDAR data exists, and propagates the analysis, as well as the derived forest height by sequentially going through the interferograms marked in "green", "yellow", "blue", "magenta", "pink", "cyan", "violet" and, finally, "grey".
The paper is organized as follows.Section 2 describes the proposed automatic mosaicking algorithm, including the nonlinear least squares representation of the fitting problem, the example of a three-scene mosaicking problem along with the simplified nonlinear least squares solution and the matrix formulation of the mosaicking algorithm generalized for multiple scenes.Section 3 demonstrates the mosaicking and validation results, as well as the associated discussions.In particular, a new mosaic map of forest height will be generated for the U.S. state of Maine through utilizing the mosaicking algorithm, compared with ground validation height maps over the Howland forest in the central Maine area and the entire state of Maine and, finally, improved with the use of a forest/non-forest map as a simple data fusion task.

Methods
In order to propagate the inverted forest height through scene overlap areas, it is necessary to explicitly define a fitting metric given any pair of overlapping forest height estimates.In particular, a nonlinear least squares problem is formulated to characterize the fitting metric.We then investigate the solution for a three-scene mosaicking problem that serves as a simplified scenario and, finally, generalize the matrix formulation for multiple overlapping scenes.

Nonlinear Least Squares Fitting Metric
To begin, a comparison is made between two sets of forest height estimates in their overlapping region.According to [15], the observed repeat-pass HV-polarized InSAR correlation magnitude, |γ HV v&t |, due to the coupled effects of volume scattering and temporal change, is related to the desired forest height estimate h v as, where S scene (unitless; 0 ≤ S scene ≤ 1) characterizes the dielectric fluctuation of the volume scatterers (perhaps due to moisture change, e.g., rainfall; a smaller S scene indicates a bigger dielectric change), while C scene (in meters; C scene > 0) represents the random motion of volume scatterers (perhaps due to wind; a smaller C scene implies a higher level of motion).Only the main lobe of the sinc function is used in Equation (1).Through inverting Equation (1), the forest height estimates can thus be considered as a function of the correlation measurements |γ HV v&t | and the fitting parameters (S scene and C scene ).Suppose there exist two sets of forest height estimates in an overlap area with the height estimates inverted as below, Here, f is the above-mentioned implicit function performing the forest height inversion, where we further omit the variable |γ HV v&t | to keep the notation concise, since the correlation magnitude is invariant in the process of data fitting.Subscripts i = 1, 2 are used to differentiate the forest height estimate, as well as the model parameters from the i-th set.
To proceed, it can be assumed that the inverted forest heights from the repeat observations, on average, are comparable to each other, and therefore, a metric is desired so that the difference between h v 1 and h v 2 can be minimized.In [15], a fitting metric comprised with two parameters was used that was comprised of the slope k and offset b.These are illustrated in Figure 2. In Figure 2, the two sets of forest height estimates are considered as the horizontal and vertical axes with the data cloud illustrated as an ellipse.The slope parameter k describes the slope of its major axis, while the offset parameter b represents the relative difference between the average forest height estimates.In particular, k and b are written as: where φ is the angle between the major axis and the horizontal axis and m 1 and m 2 are the average forest height estimates.The calculation of these parameters is obtained through a principle component analysis-based method, as described in [15].In order to have h v 1 and h v 2 match one another, a nonlinear least squares criterion is used to seek the proper model parameters S scene i and C scene i (i = 1, 2), such that the following residual error can be minimized, i.e., During a non-automated mosaicking process [15] (shown in Figure 1), for a particular overlap area, one set of forest heights is always known prior to the inversion of the other, and thus, Equation ( 6) is repeatedly used as the residual error that is minimized in order to achieve the optimal estimates of the model parameters (and, thus, forest heights) for the other InSAR scene.In other words, only the estimates from one InSAR scene are considered unknown each time the optimization step is run.However, if multiple overlapping scenes are used, it is desired to come up with the estimates simultaneously with the residual fitting error minimized globally.In order to see this, we next consider the three-scene mosaicking problem as a simplified scenario.

Three-Scene Mosaicking Problem
As shown in Figure 3, we thus apply this fitting metric (characterized by k and b) to a simple threescene mosaicking problem.There are three InSAR scenes along with a narrow validation site, where the forest heights h v 0 are predetermined and considered as ground truth data.In this example, three overlap areas can be obtained, i.e., Scene 1 with the validation site, Scene 1 with Scene 2 and Scene 1 with Scene 3. A set of three equations, f 1 , f 2 and f 3 , for these overlap areas can be summarized as: In the overlap regions (which includes the overlap of the central scene with ground validation data), two fitting parameters, i.e., k i and b i , are specified, where i ∈ {1, 2, 3}.The vector representation is defined as: where: and G is the implicit function that relates ξ to ρ.The least squares solution ρ * can be obtained by minimizing the target function: where: and " • 2 " is the Euclidean norm of a vector.Since G is a nonlinear function, we can solve this nonlinear least squares problem by utilizing the Gauss-Newton algorithm [19].Particularly, we first perform the Taylor series expansion of G at an initial point ρ 0 (that is close to the optimal point ρ * ), and keep the terms up to the first order, i.e., where J is the 6 × 6 Jacobian matrix calculated at the initial point ρ 0 and defined as: Thus, by letting ρ = ρ * and rearranging the terms in Equation ( 11), we have: where J is assumed invertible (which is usually true in practice).Because it is difficult to express the function G analytically, derivatives in Equation ( 12) are calculated numerically.If G is a linear function, the result of Equation ( 13) is exactly the desired optimal point.However, since G is nonlinear, the Gauss-Newton algorithm is an iterative numerical method, which considers the result of Equation ( 13) as a new initial point and refines the value of ρ * through another circulation of Equation ( 13).

Multi-Scene Mosaicking Problem: The Matrix Formulation for the Generalized Scenario
With the matrix form of the least squares solution for the three-scene case established, a generalized treatment of multiple connected InSAR scenes along with multiple validation sites can be formulated.To begin, out of N repeat-pass InSAR scenes that are connected with one another, M will have validation sites.The number of the connected pairs is given by E. This is shown in Figure 1, where each InSAR scene is represented as a node, each connected pair is described as an edge and only the central scene (marked in "red") has a validation site.For the state mosaic of Maine, this gives N = 37, M = 1 and E = 57.As another example, in the three-scene case of Figure 3 Further, the same vector notation as in Section 2.2 can be applied to this generalized case.In particular, Equations ( 7) and ( 9) still hold with ξ and ξ * being 2(E + M ) × 1 vectors, ρ and ρ * 2N × 1 vectors.By letting ρ = ρ * , Equation ( 11) is rewritten as: where the Jacobian matrix is 2(E + M ) rows by 2N columns.More explicitly, As long as all of the InSAR scenes are connected with one another and there is at least one validation site available, the relationship (E + M ) ≥ N holds, which satisfies the prerequisite of the Gauss-Newton algorithm [19].Therefore, in order to solve for the optimal point ρ * , the more general form of Equation ( 13) is written as [19]: Note that when J is a square matrix, Equation ( 16) reduces to Equation (13).
The manual mosaicking process in [15] sequentially propagates the information from the validation sites throughout the adjacent-scene overlap areas.Therefore, the propagation of errors will manifest itself as the "wallpapering" problem described in the Introduction to this paper.However, the nonlinear least squares solution has the benefit of estimating the desired model parameters simultaneously and minimizing the fitting error globally and, thus, resolves the "wallpapering" problem.Moreover, it is a robust system that is mathematically traceable.The accuracy and the computational complexity of this automatic mosaicking algorithm depend on the choice of the initial point, the number of iterations, the number and the quality of the InSAR scenes that are to be mosaicked together, as well as the number and the distribution of the ground validation sites.

Results
In this section, the Gauss-Newton algorithm is used to solve for the nonlinear least squares solution of the model parameters (S scene and C scene ).With these results, a new mosaic map of forest height for the U.S. state of Maine can be generated and compared with the LVIS height data and the NBCD basal area weighted (BAW) height.

Generation of the New Mosaic Map of Forest Height for Maine, U.S.
Before running the automatic mosaicking algorithm, we want to exclude one InSAR scene (i.e., "120_870" in Figure 1) and its associated pairwise connection (i.e., the directed edge pointing from "120_880" to "120_870" in Figure 1), because the temporal change effect occurring within this InSAR scene is so severe (characterized by S scene = 0.3 and C scene = 4.89 in [15]) that the global error minimization will be biased.Therefore, in this scenario, we end up with N = 36, M = 1 and E = 56.
The essential part in the implementation of Equation ( 16) is the numerical calculation of the Jacobian matrix J expressed in Equation (15).Because the elements of J are the partial derivatives of k and b with respect to S scene or C scene , numerical derivatives are calculated by allowing S scene (or C scene ) to have a small increment of 10 −6 (or 10 −5 m).The Jacobian matrix is then computed on a column-by-column basis.Given a small increment of the i-th element of ρ, the vector derivative of ξ is equivalent to the i-th column of the Jacobian matrix.The initial point ρ 0 is chosen from the model parameters determined from the manual mosaicking process [15].After ten iterations, the residual error (i.e., the Euclidean norm in Equation ( 9)) is plotted as red circles in Figure 4.It can be seen that after the third iteration, the residual error becomes very stable, which implies that the initial point (determined from the manual mosaicking process in [15]) is close to the optimal point.The effect of the initial point on convergence can be tested by assigning a uniform value (e.g., S scene = 0.65 and C scene = 13) for all of the InSAR scenes.As shown by the blue triangles in Figure 4, it is demonstrated that although the initial residual error is higher than that from the manual mosaicking process, the final result after ten iterations is similar to that when initial values were chosen from the manual mosaicking results.
By utilizing the refined model parameters (S scene and C scene ) after the tenth iteration, a new mosaic map of forest height is created for the U.S. state of Maine, as shown in Figure 5, where all of the values over water bodies have been removed by using National Land Cover Database (NLCD) 2006 [20].
A comparison of forest height estimates between the old and new mosaic maps highlights the systematic propagation of errors that occurs in the manual mosaicking process (the "wallpapering" problem).This is illustrated in Figure 6.It can be seen in the figure that the scenes that are far away from the central Maine area (where the Howland Forest is located) are more likely to have larger uncertainty in the forest height estimates (and, thus, the model parameters).Here, the absolute error in the forest height estimates between the old and new mosaic maps is illustrated and color-coded as indicated ("blue" being 0 m and "red" being 5 m).All of the values over water bodies have been removed by using NLCD2006.It can be observed that the "wallpapering" problem occurs, since the scenes that are far away from the central Maine area (where the Howland Forest is located) are more likely to have larger uncertainty in the forest height estimates.

Validation
Once constructed, a comparison can be made of the mosaic map of forest height with heights available from the LVIS sensor [18] (from the year 2009 over the Howland forest in central Maine) and NBCD BAW height [2] (from the year 2000 over the entire state of Maine).Note that it is desirable to compare this mosaic map of forest height with high spatial resolution LiDAR height data over various sites in Maine; however, this is restricted due to the LVIS data collection, which is restricted to the data strip between the Howland research forest and the Penobscot experimental forest in the state of Maine and has been done in [15].Therefore, in this work, we only compare our mosaic results with the NBCD BAW height data over the areas outside the LVIS LiDAR strip.The NBCD data consists of a 30-m resolution estimate of basal area weighted height, aboveground live dry biomass and standing carbon stock for the conterminous United States in 2000.The dataset was derived by utilizing the empirical modeling approach, which combines USDA Forest Service Forest Inventory and Analysis (FIA) data with 2000 SRTM InSAR data and optical remote sensing data acquired from the Landsat ETM+ sensor [2].

Comparison with LVIS and NBCD BAW Heights over Howland Forest
Although the Maine mosaic of forest height (i.e., Figure 5) was determined with the use of the LVIS height at the ground validation site, it is useful to compare the inverted height estimates from the final mosaic with the LVIS data.As shown in Figure 7, it can be seen that the ALOS InSAR correlation magnitude-inverted heights correspond visually well with the LVIS heights.By contrast, the NBCD BAW height seems to have a much shorter dynamic range, although it does indicate similar features as those derived from LVIS and ALOS data.Note that the underestimation (termed as "inherent bias" in this work) of the NBCD BAW height can be observed in comparison to the LVIS height data, especially over the urban area that is marked by a "red" window in Figure 7d.A quantitative comparison of results is illustrated in Figure 8.Each point in Figure 8 represents a forest area of 400 m × 800 m derived from multi-pixel averaging.Figure 8a shows the consistency between the ALOS InSAR correlation magnitude-inverted heights and the LVIS heights.The root mean square error (RMSE) of this fit is 3.8 m, and the statistical R value is 0.48. Figure 8b compares the NBCD BAW and LVIS-derived heights.In this part of the figure, it can be seen that the NBCD BAW height has a smaller dynamic range, as well as some inherent bias compared to the LVIS heights.This is likely due to the fact that the NBCD is basal area weighted and less sensitive to the height of dominant trees within a resolution element, as is the case for the LVIS-derived heights.The presence of unusually high temporal decorrelation (e.g., in agriculture and water-covered regions) is known to bias forest heights derived from the ALOS InSAR correlation measurements.This can be seen in Figure 8a at the lower end of the LVIS-derived height range, where water bodies, farmlands and urban activities, not thoroughly removed from the ALOS imagery, yield forest height estimates significantly larger than those observed by LVIS.Regions such as this can be detected because of the estimates of unrealistically large trees (40 m and larger) and/or removed in the larger mosaic using a forest/non-forest classification map, which will be discussed later.

Comparison with NBCD Heights over the Entire State of Maine
The mosaic map of the NBCD BAW height is illustrated in Figure 9 for the entire state of Maine.Compared to Figure 5, it can be seen that the mosaic of NBCD BAW height has a much shorter dynamic range.The quantitative comparison result between the mosaic of forest height inverted from ALOS InSAR correlation magnitude (i.e., Figure 5) and the mosaic of the NBCD BAW height (i.e., Figure 9) is shown in Figure 10.Each point in Figure 10 corresponds to a forest area of 500 m × 500 m through multi-pixel averaging.The inherent bias of the NBCD BAW height that appears in Figure 8b can also be seen here for the whole mosaic map.Overestimation of forest height from the ALOS InSAR correlation magnitude-inverted height, due to the temporal decorrelation in urban areas, farmlands etc., is evident in the upper right-hand corner of Figure 5.This part of the state of Maine consists of scattered farmlands and shows up as the colors "orange" and "red" in the imagery, indicating heights of 35 m and taller, which is much larger than average tree heights in the region.While this may be a useful tool for detecting change, even within a forest (e.g., selective logging), here, it is considered a primary source of error, which can be improved by combining the mosaic with a land cover database that differentiates forested and non-forested regions.Typical forest/non-forest maps have already been derived from the ALOS SAR backscatter power, as demonstrated in [7,11].However, as noticed in this work, another resource for such a classification can be the NBCD BAW height mosaic, where any non-forest region is identified with the use of a flag value.The refined mosaic map is illustrated in Figure 11.Comparing with Figure 5, it can be seen that the overestimated height values over the non-forest regions (shown as "orange" and "red" spots in Figure 5) have been removed in the updated forest height mosaic.The refined mosaic map after removing the height estimates over non-forest regions by using the NBCD mosaic.This mosaic is also color-coded as indicated ("blue" being 0 m and "red" being 45 m).All of the values over water bodies have been removed by using NLCD2006.

Discussion
In this work, there are several practical concerns that should be considered with the proposed improvements that are related to the implementation of this forest height inversion approach and its automated mosaicking process.
First, temporal decorrelation (e.g., harvesting over farmlands, urban activities) is noticeable and embodied as overestimated height values over the non-forest regions in the mosaic map.One practical approach is to remove the non-forest regions through the use of a forest/non-forest map.For example, in this work, the NBCD mosaic map has been utilized to serve as a forest/non-forest map.However, forest/non-forest maps derived from SAR backscatter power [7,11] can alternately be used so that SAR/InSAR observations from the same spaceborne mission are fully exploited (i.e., InSAR correlation magnitude data are used to generate a forest height mosaic, while SAR backscatter power is used to create a forest/non-forest map).This would alleviate the need for external maps.Further, a forest/non-forest classification map can be applied to the InSAR coherence map prior to the forest height inversion, instead of being a post-processing step to the mosaic results, as in this paper, so that the accuracy of the forest height inversion can thus be improved by precluding non-forest regions being used in overlap regions.Note that the forested plots that are affected by selective logging and/or forest degradation cannot be removed by using the forest/non-forest map and will also embody themselves as overestimated "large" forest heights, which could be an interesting result all by itself and useful for monitoring the global forest change.
Second, because of the repeatable nature of SAR data collections, there are often many scenes available over the same area, but separated in time by weeks, if not months.However, due to the unreliable nature of the temporal decorrelation effects, only a few of them are suitable for the use of forest height inversion.Compared to a stable weather condition, a windy and/or rainy day will decrease the observed InSAR correlation magnitude by a great amount.Although the data with smaller correlation magnitude still have the vegetation structural and temporal change information that could be utilized for forest height inversion, this bit of information is often masked by correlation sampling noise [21], making the inversion much noisier and less robust.Furthermore, if the weather condition changes non-uniformly, such as a regional rainfall, the temporal change effects may vary across each InSAR scene, so that the model parameters cannot be assumed constant over the whole scene any more.In this paper, and in previous work [15,22], through a careful selection of ALOS InSAR scenes over the same study area, only one or two out of the dozens of available scenes are best suited for forest height inversion.It is recommended and desired to have more reliable spaceborne repeat-pass InSAR data with moderate (less than a month; 12 days for NISAR [23] and 14 days for ALOS-2 [24]) or large (on the order of months; 46 days for ALOS) temporal baselines, so that the best InSAR scene(s) can be selected and utilized to generate a reliable forest height mosaic.

Conclusions
In this paper, an automatic mosaicking algorithm for creating large-scale mosaics of forest height using spaceborne repeat-pass cross-pol InSAR correlation magnitude data is demonstrated.In order to invert forest height from repeat-pass InSAR data, two model parameters (S scene and C scene ) are estimated for each InSAR scene.This method improves on a manual mosaicking approach [15] used previously, that propagates information and error from the validation sites throughout adjacent scenes.This manual mosaicking approach creates an effect of increasing errors as the distance of estimates from the validation site increases, an effect termed as the "wallpapering" problem.
The automated method derived in this paper has the benefit of equalizing errors across the mosaic, simultaneously, and minimizing them globally.This is achieved by formulating two fitting parameters (slope k and offset b), which describe the differences between the height estimates from two overlapping InSAR scenes, and minimizing the differences using a nonlinear least squares fitting routine.In order to determine the global minimum of this nonlinear least squares problem with multiple InSAR scenes, a Gauss-Newton successive minimization algorithm is used.The proposed mosaicking algorithm was applied to 36 ALOS InSAR scenes, with relatively high correlation magnitude and uniform model parameters across each scene, and a new mosaic map of forest height was generated for the U.S. state of Maine.It was noticed that the convergence of the model parameters was achieved after the third iteration of the Gauss-Newton algorithm.
This forest height mosaic was compared with both LVIS and NBCD BAW heights over the Howland Forest in the central Maine area and also compared with the NBCD mosaic over the entire state of Maine.The validation results show that the ALOS InSAR correlation magnitude-inverted heights correspond with the LVIS and NBCD BAW heights; however, the NBCD BAW heights have a much shorter dynamic range and are mostly limited to heights below 15 m.Finally, the ALOS InSAR correlation-derived mosaic map of forest height was refined by utilizing the NBCD mosaic as a forest/non-forest map.
The potential limitations about this forest height inversion approach and its mosaicking method involve the temporal decorrelation over non-forest regions, as well as the forest disturbance, such as selective logging.As discussed in Section 3.3, this limitation could be partially avoided by incorporating a forest/non-forest classification map to mask out the non-forest regions, while the forest disturbance cannot be removed (however, this shows a way for forest change detection, an interesting result all by itself).Another limitation would be the unreliable nature of the temporal change effects in repeat-pass InSAR observations.This error could be reduced through the collection of more spaceborne repeat-pass InSAR observations made with a shorter revisit period (less than the ALOS-1 46-day period).With shorter revisit periods, further improvement could be achieved through the incorporation of a volumetric component in the InSAR model.This would require a non-zero baseline and would nominally be balanced against the temporal decorrelation signature used here.
By noting and managing these limitations properly, the techniques proposed in this paper in combination with its companion paper [15] are an efficient tool for creating a state-and/or continental-scale mosaic of forest height from spaceborne repeat-pass cross-pol InSAR correlation magnitude.Importantly, this method demonstrates an approach for using archival repeat-pass InSAR observations (e.g., JERS-1, ALOS, ALOS-2, NISAR) to map vegetation height over large regions, potentially at a continental scale, and, to our knowledge, is the first method that successfully utilizes spaceborne repeat-pass InSAR data to create large-scale forest height mosaics in the InSAR vegetation community.The algorithms described here can be implemented with reasonable computational cost and, thus, serve as an observational prototype for repeat-pass InSAR missions, like ALOS-2 and NISAR.The present approach also serves as an alternative and complementary tool for other PolInSAR inversion techniques when single-pass InSAR and/or full-polarization data may not be available.

2 Figure 2 .
Figure 2. Geometric illustration of the comparison between two sets of forest height estimates.The data cloud is considered as an ellipse, with the angle between the major axis and the horizontal axis denoted by φ and their average heights denoted as m 1 and m 2 .

Figure 3 .
Figure 3. Illustration of mosaicking three InSAR scenes.The "blue" bar is a ground validation site with the heights predetermined.In this example, there are three InSAR scenes (i.e., three pairs of model parameters S scene and C scene need to be determined) along with three overlap areas (i.e., three pairs of fitting parameters k and b can be computed).
, N = 3, M = 1 and E = 2.For the generalized scenario, there are (E + M ) overlap areas (and, thus, 2(E + M ) fitting parameters k and b) in total along with N InSAR scenes (and, thus, 2N model parameters S scene and C scene ).

Figure 4 .Figure 5 .
Figure 4. Residual error at each iteration of running the automatic mosaicking algorithm.The blue triangles indicate the results by using the average values of the model parameters (S scene = 0.65 and C scene = 13) for all of the InSAR scenes as the initial point, while the red circles show the results by considering the model parameters determined from the manual mosaicking process[15] as the initial point.

Figure 6 .
Figure 6.Illustration of the residual error between the automated and manual mosaicking approaches.Here, the absolute error in the forest height estimates between the old and new mosaic maps is illustrated and color-coded as indicated ("blue" being 0 m and "red" being 5 m).All of the values over water bodies have been removed by using NLCD2006.It can be observed that the "wallpapering" problem occurs, since the scenes that are far away from the central Maine area (where the Howland Forest is located) are more likely to have larger uncertainty in the forest height estimates.

Figure 7 .
Figure 7.The (a) optical image from Google Earth is compared with (b) the LVIS height, (c) the ALOS InSAR correlation magnitude-inverted forest height and (d) the National Biomass and Carbon Dataset (NBCD) basal area weighted (BAW) height over the Howland research forest in central Maine.The values over water bodies are removed with the use of National Land Cover Database (NLCD) 2006.The color maps are coded with the same color scale ("blue" being surfaces, "red" being 45 m) and a spatial resolution of 50 m × 50 m.The inherent bias of the NBCD BAW data, discussed in the text, is highlighted by a "red" rectangular window over the urban area.

Figure 8 .
Figure 8. Quantitative comparison results between various height metrics.(a) The comparison between the ALOS InSAR correlation magnitude-inverted height (Figure 7c) and the LVIS height (Figure 7b) with RMSE of 3.8 m and R value of 0.48; (b) the comparison between the NBCD BAW height (Figure 7d) and the LVIS height (Figure 7b) with RMSE of 5.6 m and R value of 0.3.Each point corresponds to a forest area of 400 m × 800 m through multi-pixel averaging.The data points pertaining to the inherent bias of the NBCD BAW height are indicated by a dashed circle.

Figure 9 .
Figure9.The mosaic map of the NBCD BAW height for the state of Maine, U.S.This mosaic is also color-coded as indicated ("blue" being 0 m and "red" being 45 m).All of the values over water bodies have been removed by using NLCD2006.

Figure 10 .
Figure 10.Quantitative comparison result between the mosaic of forest height inverted from ALOS InSAR correlation magnitude (i.e., Figure5) and the mosaic of the NBCD BAW height (i.e., Figure9) for the entire state of Maine, U.S.Each point corresponds to a 500 m × 500 m forest area through multi-pixel averaging.The data points that are affected by the inherent bias of the NBCD BAW height and by the temporal decorrelation of the ALOS InSAR data (e.g., farmlands and urban area) are indicated by dashed circles, respectively.

Figure 11 .
Figure 11.The refined mosaic map after removing the height estimates over non-forest regions by using the NBCD mosaic.This mosaic is also color-coded as indicated ("blue" being 0 m and "red" being 45 m).All of the values over water bodies have been removed by using NLCD2006.