An Algorithm Measuring Urban Building Heights by Combining the PS-InSAR Technique and Two-Stage Programming Approach Framework

Urban building heights offer critical information in studying urbanization. An approach with adequate accuracy is of significant interest. Spaceborne interferometric synthetic aperture radar (InSAR) techniques and multibaseline (MB) InSAR observations can be an alternative to fill the interest. Since the MB InSAR observations are achieved through multitemporal SAR observations with different baselines, temporal decorrelation is unavoidable, affecting most existing MB-based phase unwrapping (PU) approaches. Also, the quality of the derived heights depends on a PU approach. Thus, by combing the permanent scatterers (PS) InSAR technique and the two-stage programming approach (TSPA) framework, we proposed an algorithm to estimate urban building heights with the multitemporal SAR observations, calling it the PS-TSPA algorithm. It takes advantage of the global phase gradient information used in the single baseline PU approaches, not considering the phase continuity assumption. The algorithm should be robust to the temporal decorrelation or noise effect. With the time series of 17 TerraSAR-X data, we estimated building heights in the southern urban area of Chengdu, China. The derived building heights were satisfactory overall and assessed by in situ measurements. In comparison, the Tomo-PS method estimated the heights in the same area and with the same SAR data. Quantitative analyses indicate that the proposed algorithm outperforms the Tomo-PS method. Finally, with routine SAR observations and SAR data downloadable through the Internet, the PS-TSPA algorithm can be inexpensive in estimating urban building heights with reasonable accuracy.


I. INTRODUCTION
A BUILDING height in an urban area is the most intuitive information about human concentration, population density, and commerce. It is vital for urban research, including the population [1], energy consumption [2], and urban expansion [3], [4]. With the successful launching of spaceborne synthetic aperture radar (SAR) and widely available multiobservation SAR datasets through Internet downloading, spaceborne remote sensing data and techniques contribute to and even revolutionize urban mapping research (e.g., [5]). The SAR data become indispensable for their characteristics of multiple observations and routine operations [6], [7], [8]. SAR, an active remote sensing sensor, is weatherindependent. It is particularly advantageous for studying a city (e.g., Chengdu, China) with constant and pervasive cloud/fog cover. Various SAR-based techniques, including the intensitybased, 3-D tomographic SAR (TomoSAR) [9], [10], [11] and interferometric SAR (InSAR) methods [12], [13], [14], [15], estimate building heights in urban areas. The intensity-based method generally requires less SAR data, but its success is affected by the building's geometric shapes and materials and proper SAR imaging geometry. The shapes and materials are countless, but the imaging geometry is usually fixed. Thus, obtaining accurate building heights in urban areas based on radar backscatter intensity is quite challenging [8], [16].
The TomoSAR methods obtain 3-D information by forming a synthetic aperture in the elevation direction with stacks of repeat pass SAR observations or dataset, whose core is to estimate the number of permanent scatterers (PS) in overlapping pixels [11]. The methods can usually be divided into the spectrum estimation [17], [18] and compressed sensing (CS) [19] methods. The spectrum estimation methods can be noise-prone and generally low in elevation resolution and accuracy (e.g., [11]). The CS algorithm considers the signal existing sparsely, seeking high elevation resolution. The mathematical principle can be expressed as y = Aγ, where y is a column vector with N elements or observations. A is the N × M sensing matrix. γ is another column vector with M elements, or it is a sparse signal to be solved, i.e., tomography, pixel by pixel. Since N is much less than M, y = Aγ is ill-posed, and one needs to solve a nonlinear optimization problem to obtain an approximation for the signal. The process is generally complicated and computationally intensive. Therefore, designing a CS-based TomoSAR algorithm considering the accuracy and computational efficiency [20] is challenging.
An InSAR method is valuable for surface elevation estimation [21], [22], [23]. Its data acquisition requirement, as compared to the TomoSAR data acquisition, is much easier met. In addition, with PS dispersed in a study area, influences of the temporal and geometrical decorrelations and atmospheric phase screen can be reduced [24], [25]. Furthermore, elevation refinement with an external digital elevation model (DEM) for PS can estimate a target or building height. Nevertheless, the success of refining the PS-related methods depends mainly on a phase unwrapping (PU) algorithm to unwrap the phase of each permanent scatterer in the urban area, which is another significant challenge.
According to the number of interferograms used, the PU algorithms can be divided into the single-baseline (SB) or multibaseline (MB) ones. The former assumes the phase continuity (i.e., the Itoh condition) to obtain the absolute phase [26]. Unfortunately, the assumption is usually not satisfied in building height estimation due to the rapid phase change or significant phase gradient from the bottom to the top of a building or skyscraper. The inadequate SAR sampling rate that does not satisfy the Nyquist-Shannon sampling theorem is attributed to be the primary cause. The MB PU techniques can relieve the Itoh condition requirement with perpendicular baseline (hereafter referred to as the baseline) diversity. The MB PU techniques can potentially achieve satisfactory PU results in an area with buildings/skyscrapers [27]. Nevertheless, most MB PU methods (e.g., [28], [29]) are noise-prone, which makes them uneasy about being effective and routine in use.
Although Yu and Lan [30] developed a two-stage programming approach (TSPA) framework, reducing the noise's effect, the framework was not intended for estimating a building height because a regular-shaped grid with rectangular cells is used in PU. On the other hand, many PS exist in an urban area, and they are usually dispersed. In building height estimation, one has to use the scatterers that form a network with an irregular grid. The revision of the regular-shaped grid in [30] is in order. Thus, to take advantage of ample PS in an urban area, we combine the InSAR permanent scatterer model with the TSPA framework, calling the combination a PS-TSPA PU algorithm. Procedurally, the algorithm consists of two primary steps. First, a spatial network (e.g., the Delaunay network) approach establishes the adjacency relationship of neighboring PS. The Chinese remainder theorem (CRT)-based method (e.g., [30]) is used to estimate the gradient of ambiguity numbers of the neighboring PS to mitigate the effect of not meeting the Itoh condition, recognized as one primary deficiency in the SB PU approaches. Second, the ambiguity number in the PS-TSPA PU algorithm is calculated through global optimization. Thus, the algorithm's noise robustness is enhanced, reducing the noiseprone deficiency of existing CRT-based MB PU algorithms. With unwrapped phases, we then estimate building heights. Since we focus on the permanent scatterer's PU problem in the workflow, the analysis is simplified and is convenient for practical applications. To assess the algorithm's effectiveness, we use 17 TerraSAR-X data, deriving building heights in urban areas of southern Chengdu, China. As shown in the study, the algorithm outputs building heights satisfactorily and effectively.

A. Basic Principle of the InSAR Analysis
A permanent scatterer and InSAR phase unwrapping: The scatterer within a dataset of multiple registered SAR images is usually selected if the pixel has a small amplitude dispersion index (ADI) and a high coherence coefficient. It is characterized by high and stable radar backscatter. A building or an exposed rock is a good candidate for the permanent scatterer. Achieving satisfactory elevation or surface deformation estimation requires multiple and dispersed PS.
The InSAR PU of a selected pixel, s, is where ϕ(s) is the wrapped phase, ψ(s)is the absolute phase, k(s) is an ambiguous integer, and ψ(s) and k(s) are unknown. Thus, (1) is ill-posed, and additional information is needed to obtain a unique solution.
An SB PU algorithm for PS: It generally consists of two steps [31]. First, the Delaunay or other networks link the neighboring PS, forming elementary closing loops. The phase gradient of two adjacent PS, u and v, is estimated along the arc linking both scatterers. Second, the absolute phase gradient of u and v is solved by making it as close as possible to the estimated gradient using an L p -norm optimization model. Therefore, the estimated phase gradient generally determines the accuracy of PU results.
Mathematically, the SB PU method utilizes the phase continuity assumption to estimate the gradient of u and v, i.e., the absolute gradient value of neighboring PS is less than π [26] or whereΔψ(u, v) is the estimated phase gradient. After obtaining the gradient, the L p -norm optimization model minimizes the difference between the absolute and estimated phase gradients. The minimization process [32] is where n represents the number of pairs of neighboring PS, i.e., the number of the network arcs, ω is the weight, and t is the auxiliary variable. As previously stated, if conditions of (2) are not met, the PU uncertainty or error occurs. An MB PU algorithm: When two flattened interferograms are considered in phase unwrapping, it is called a double baselines (DB) PU algorithm, or is the ith baseline length of the ith interferogram, and B 1 and B 2 differ in length. ϕ i (s)is the ith wrapped phase and k i (s) is the ith ambiguous integer. When m (m = 3) baselines are considered, one extends (4) as follows: Similarly, one can write the PU algorithm when m > 3. The MB PU methods simultaneously solve integers k i (s), (i = 1, 2, ..., m) (e.g., [27]). Then, the absolute phase can be seriously affected by noise since the global gradient phase information is not used. However, the information is typically used in an SB PU algorithm to reduce the noise effect significantly.

B. Proposed PS-TSPA PU Algorithm and a Building Height Estimation
The overarching objective is to develop a new MB PU algorithm for PS within a study area. The algorithm takes advantage of the global phase gradient information but does not consider the phase continuity assumption. Thus, it is robust to the noise effect and obtains satisfactory phase unwrapping results and building height estimation with the absolute phases of PS.
The PS-TSPA PU algorithm: Considering previous discussions of the PS, the SB and MB PU approaches, and the TSPA, we express the proposed algorithm with m = 3 as an example.
(The algorithm with m > 3 can be analogically obtained.) Mathematically, the algorithm for two adjacent PS is Subtracting the two equations in (6), one obtains the following: where is the gradient of ambiguity numbers of the two scatterers. Since k i (u) and k i (v) are integers, Δk i (u, v) is an integer too and can be estimated [30].
With the estimated gradients of the two adjacent PS, the PU results of the m interferograms can be separately derived with the L p -norm optimization model, like in an SB PU method. The model for the ith interferogram is where n represents the number of pairs of neighboring PS, i.e., the number of the network arcs.
In (3) and (8), different p values correspond to variable PU methods. When p = 1 and 2, the related methods are the minimum-cost flow [33] and the least square one [34], respectively. When p = 0, (3) and (8) become NP-hard problems, where the NP stands for nondeterministic polynomial. Then, (3) or (8) cannot be solved in polynomial time. Alternatively, in the SB PU domain, the branch-cut algorithm [35] can obtain an approximate solution of the L 0 -norm. In the MB PU domain, (8) can be approximatively solved by the outlier-detection solution method [36], [37], which iteratively uses the L 2 -norm method to remove the constraints in (8). When the ambiguity number of each permanent scatterer is solved, its absolute phase is obtained by (1).

Estimating a building height in an urban area:
The radar wave usually hits the surface like a building rooftop first in the earth observation with a spaceborne sensor. The rooftops usually have short dividing or side walls, air conditioning units, and metallic vents or ducts. They interact with the rooftop surface, forming corner reflectors. As the rooftops do not vary much over time, the reflectors possess stable radar backscatter and show high coherence in an interferogram. The reflectors are then considered PS, forming the digital surface model (DSM) in the InSAR analysis for an urban study. With the help of accurate temporal PU results obtained by the proposed PS-TSPA PU algorithm, building heights are obtained by calculating the PS heights after removing the non-height-related phase errors caused by variable atmosphere conditions and system noise.
To estimate building heights in an urban area, we remove the terrain phase (ψ topo (s)) and atmospheric phase (ψ atmo (s)). Then, the residual phases contain the building height (ψ BH (s)) information. Thus, for the ith interferogram, the building height (h) related to one permanent scatterer (s) is is the phase of the linear deformation along the line-of-sight direction, v(s) is the surface deformation rate (m/day), and T(i) is the number of days between two acquisition dates of SAR data forming the interferogram. With the leastsquared fitting to (9), h is obtained. Fig. 1, the flowchart of the proposed PS-TSPA PU algorithm to estimate a building height, consists of seven primary steps.
1) Selecting one SAR image as the common reference (i.e., the primary image) and geo-registering the rest images to it, one by one. 2) Forming interferograms. To reduce the impact of nonlinear deformation and atmospheric phase on the unwrapping phase results, we chronically create interferometric pair. 3) Selecting PS and forming a permanent scatterer network. 4) Removing the topographic phase with an external DEM. 5) Deriving the time series of absolute phases of each permanent scatterer. 6) Filtering nonheight related phases caused by atmospheric variations. 7) Fitting the residual phase of each permanent scatterer and estimating the building height. We call the PS-TSPA PU algorithm and height estimation collectively the PS-TSPA algorithm onward for short.

A. Chengdu and its Southern Urban Area
Chengdu City is the capital of Sichuan Province, located in southwestern China. The city has an estimated population of ∼10 million in core urban areas in 2023, with an area of 14 380 km 2 . Over one-half of the Fortune 500 companies have offices here. The city, the most innovative in China and the world, is a leading hub city in China in research and development (R&D) and commerce (https://worldpopulationreview. com/world-cities/chengdu-population).
Chengdu, located in the western part of the Sichuan basin, has a subtropical monsoon climate with annual precipitation of up to 1000 mm. The precipitation is mainly rainfall, falling the most in summer and sufficient in winter. With the basin topography, the air is moist year-round. The cloud/fog cover is annually pervasive with spotty sunny days. Any SAR-based algorithms to estimate building heights in the city are preferred. Moreover, many high-rise office and residential buildings exist, mainly in the southern urban area of Chengdu, where the phase continuity assumption is hardly met. It is a significant challenge to the SAR-observation-based technique. Fig. 2(a) is a GoogleEarth image acquired on 26 November 2016, showing the study of about 5.5 km × 7.6 km (The image was flipped east-west geographically. As discussed next, the SAR images were flipped east-west since they were acquired along descending polar orbits with the SAR sensor looking right and down toward the ground. Thus, both types of images are oriented roughly similarly). The east or left side area has many regular/irregular-shaped tall buildings with high densities. Building height measurement from the ground is generally challenging due to lacking adequate open space around each building. Four buildings (A-D) are identified as validation targets for the study (Section III-C). Buildings on the west or right side are generally not tall, with a sparse building distribution. For instance, the maintenance facilities of the Taipingsi Airport and warehouses and factories around are low in height. The airport, built over 80 years ago, was small and is out of commission. Instead, the Shuangliu International Airport of Chengdu is in operation and off the image southwest.

B. TerraSAR-X Data and Interferograms
Seventeen TerraSAR-X single look complex data acquired along descending orbits between 14 March 2016 and 13 August 2017 are available and are studied. The spatial resolutions in the azimuth and slant range directions are 1.98 and 1.36 m, respectively. Then, a subset operation was performed on the seventeen images, creating 17 subimages. Each subimage has 3856 rows (azimuth) by 4010 columns (range). Fig. 2(b) shows the average intensity of the 17 images in greyscale. Since SAR uses a side-looking geometry to image the ground, geometric distortions (e.g., overlay, foreshortening, and shadows) occur. The distortions pose a significant challenge in building height estimation using SAR observations and InSAR techniques.
Considering the procedure of selecting a reference image in the PS-InSAR technique [25], [38] and among the 17 subimages, we chose the SAR image acquired on November 11, 2016 as a reference or a common primary image. The rest 16 SAR images or secondary ones were registered to it, one by one. Then, the interferometric pair was chronically made, i.e.  Table I). The baseline of each pair is also given in the table. The DEM with a spatial resolution of 12.5 by 12.5 m was used to remove the topography.

C. Four Building Heights
Standing on the ground with a laser rangefinder, Rasger T1800BE, we measured building heights. The rangefinder has a measurement error of ±1 m and a working distance from 6 to 1500 m. The heights of four tall buildings labeled A--D near Tianfu Avenue (see Fig. 2 open space around each high building is barely available. A highway/road may be an "open space" in most locations. Due to a high traffic volume, the in situ data collection is challenging and unsafe.)

A. Building Height Estimation
According to the workflow in Section II-B and with the ADI value (0.56) and average coherence value (0.5) as thresholds, 339 063 PS were identified within the time-series SAR dataset. The results are shown in Fig. 3(a). The scatterers were dispersed, with a slightly higher density on the left-and-up quad. The Delaunay triangulation algorithm connects them, forming a network of triangles with variable shapes and sizes. Fig. 3(b) is the network. The arcs, especially short ones, are clustered in the left-and-up region. The phase gradient of two neighboring PS is calculated along the connecting arc.
Then, the proposed PS-TSPA algorithm was used to obtain the absolute phase of PS simultaneously, and the non-height-related phases were removed through filtering [25]. In reference to the DEM, the DSM information for each permanent scatterer, i.e., urban target, was obtained by the least-squared fitting with the time series unwrapped phase. The distribution of urban target heights derived from the PS-TSPA algorithm is shown in Fig. 4(a). The PS in red or blue represent tall or short targets. As shown in the figure, most tall buildings are distributed on the left side and from the middle to the bottom. Based on Fig. 2(a) and on-site observation, the area, from the Tianfu Overpass with multilayer highways to the bottom, is clustered with many commercial centers with high-rise multipurpose buildings. The rest area of Fig. 4(a) consists of yellow, yellowish-green, and blue dots, indicating low height values. For instance, from the overpass and upward, there primarily are the South Railway Station of Chengdu, roads, residential areas, and low office buildings. The lowest heights are near the Taipingsi Airport, the third ring road, and the expressway linked to the Shuangliu International Airport. In short, the derived building heights are encouraging and satisfactory overall.
Intuitively, the PS-InSAR technique can be viewed as a particular case of the TomoSAR technique if a single permanent scatterer exists per pixel. Combining the TomoSAR and  [39] developed the Tomo-PS InSAR method, mapping building heights in Hong Kong, China. Using the Tomo-PS InSAR method, we also mapped building heights with the same 17 TerraSAR-X data in the study area, showing the results in Fig. 4(b). Although Fig. 4(b) and (a) may be consistent with each other, height differences exist, as shown in Fig. 5(a), the [ Fig. 4(a)-(b)] image. The differences occur mainly on the left side, where a wide range of urban targets (e.g., the railway station, road surfaces, overpasses, and low and high buildings) with different heights exist. The boxplot of Fig. 5(a) shows the distribution with a lower quartile of 1.9 m, a median of 9.8 m, and an upper quartile of 17.8 m. The mean is 14.1 m, and one standard deviation is 29.9 m. Thus, the difference should be quantitatively significant.

B. Validation
We assess both methods with four in situ building heights, indicating which one may be better in this study. Fig. 6(a) shows an office building (AOB) labeled A in Fig. 2. With ample public open spaces and narrow vegetation strips around the building, the environmental impact on the side-looking SAR image geometry is minimal. In addition, the building can be easily identified on the SAR images due to its relatively simple surroundings. Fig. 6(b) shows the heights of PS related to the building, showing different heights. Boxplots of the heights are plotted in Fig. 6(d). The minimum, lower quartile, median, upper quartile, maximum, mean, and one standard deviation are 16.0, 34.7, 49.1, 63.5, 86.6, 48.3, and 15.8, respectively (see Table III). For comparison, Fig. 6(c) shows the permanent scatterer heights derived by the Tomo-PS method, and Fig. 6(d) is the boxplot of the heights. The descriptive statistics are also tabulated in Table III. The in situ building height is 82.8 m (see Table II). If the maximum value represents the building height, the building is 86.6 m tall, derived by the PS-TSPA algorithm. It becomes 80.2 m high by the Tomo-PS method. Although both heights should be comparable, the PS-TSPA algorithm should offer a closer building height measurement overall than the Tomo-PS method due to a smaller measured range and one standard deviation. A small variation should indicate a better-converged measurement. Fig. 7(a) is the Tianhe International Community (TIC), B in Fig. 2. According to the on-site observation, the community is residential, consisting of a tall building (pointed by a white arrow) and a short one. The tall building is 84.8 m high (see Table II), but the height of the short building cannot be measured from the ground due to a lack of publically accessible open spaces or blocking by the tall building. Fig. 7(b) displays the permanent scatterer heights of both buildings after the PS-TSPA algorithm, and Fig. 7(d) is the boxplot of the heights above the dotted red line, indicating the height of the tall building.  The descriptive statistics are given in Table III. Comparatively, Fig. 7(c) shows the heights of PS derived by the Tomo-PS method, and the boxplot of the heights (above the dotted line) is shown in Fig. 7(d). The descriptive statistics are also given in   Table III. The range and one standard deviation by the PS-TSPA algorithm are still smaller than those by the Tomo-PS method, indicating a better height measurement. In addition, the highest value for the short building in Fig. 7(c) is 138.0 m after the Tomo-PS method, an obvious concern. The highest value of the short building after the PS-TSPA is 65.9 m, about 20 m shorter than the tall building.
In Fig. 8, the Chengdu Archive (CA) building (C in Fig. 2) is on the right of the red dotted line. The buildings on the left are gated residential without public access. Per the on-site observation, the archive building is taller than the residential buildings. Fig. 8(b) and (c) are the permanent scatterer heights of the area derived from the PS-TSPA and Tomo-PS methods, respectively. The scatterers from the archive and residential buildings are mixed, which can affect the height estimation. Boxplots of the heights of the archive building after both methods are shown in Fig. 8(d), and the corresponding descriptive statistics are given in Table III. Again, the range and one standard deviation by the PS-TSPA algorithm are smaller than those by the Tomo-PS method. Then, heights derived by the PS-TSPA algorithm vary less than the Tomo-PS method. The archive building is 71.4 m tall (see Table II). It is 72.9 m tall after the PS-TSPA algorithm and becomes 160.9 m after the Tomo-PS method.
The Innovation Time Center (ITC) (D in Fig. 2) area consists of four high-rise multipurpose buildings [see Fig. 9(a)]. The building in the middle, pointed by a white arrow, is of interest due to the impact of the buildings on each side. With the polar-orbit and side-looking SAR observations, the shadowing, overlay, and foreshortening should be the most severe for the building of interest, possessing a significant challenge to estimate the building height with the SAR observation and InSAR technique. The heights derived from the PS-TSPA and Tomo-PS methods are shown in Fig. 9(b) and (c), respectively. Due to the impact of the nearby buildings, we could only outline a general area representing the building of interest and extract PS representing the building heights. Fig. 9(d) shows boxplots of the heights within the white rectangle. In the plots, there are more outliers for the Tomo-PS measurements, suggesting more height jumps caused by not meeting the phase continuity assumption. The descriptive statistics are summarized in Table III. The range and one standard deviation by the PS-TSPA algorithm are still smaller than those by the Tomo-PS method. The building of interest is 134.8 m high (see Table II

A. Applicability to a Large-Scale Infrastructure
The South Railway Station of Chengdu, representing a large-scale infrastructure, is within the study area, as shown in Fig. 10(a). The image was flipped east-west. The platform and tracks are roughly horizontally oriented. Fig. 10(b) shows the heights of PS after the PS-TSPA algorithm. The spatial distribution of measured heights is generally consistent with on-site observation. The left side is higher, reaching its highest in the center of the station building and decreasing rapidly to the right side. The height boundary, pointed by a white arrow, between the station roof and railway tracks is evident when compared to the optical image [e.g., Fig. 10  Similarly, the heights of PS of the railway station area were derived from the Tomo-PS method, and the results are shown in Fig. 10(c). Most height values were negative. For example, the station building rooftop is -7.1 m, which is unacceptable. The height boundary between the rooftop and railway tracks is not apparent. Then, we extract the height of the same permanent scatterer from the track and the height of the same permanent scatterer from the rooftop. The height difference is 1.5 m. The result should not be intuitive. Boxplots of derived heights are shown in Fig. 10(d). The height ranges from -27.3 to 49.1 m after the PS-TSPA algorithm, but the range is between -39.5 and 7.2 m after the Tomo-PS method. Incidentally, 12% of PS are negative in height values after the PS-TSPA algorithm, but the percentage is 98% after the Tomo-PS method. Thus, the PS-TSPA algorithm should be preferred in estimating heights for large-scale infrastructure facilities.

B. Possible Deficiency in the Height Estimation
Lack of identifiable PS from a building rooftop: SAR observes a ground target based on its backscatter. For a mono-state and side-looking SAR sensor, the backscatter is none or weak when the surface is smooth or slightly rough because the incidence radar signal energy is mostly specularly reflected. Without adequate backscatter from the rooftop, PS representing the rooftop heights lack. Then, the height estimation becomes unavailable. To support the argument, we example the Teda Times Center (E in Fig. 2) area that consists of several high-rise multipurpose buildings [see Fig. 11(a)]. A white arrow points to the building of interest. Fig. 11(b) and (c) are PS representing building heights in the area. However, the scatterers related to the building of interest may lack, or the scatterer delineation is uncertain. Thus, the height estimation is affected. Caution should be exercised.
The irregular network formed by PS: As shown in Fig. 3(b), long arcs exist when the ground target has a large and flat extent [e.g., the Taipingsi Airport and nearby roads, outlined by a rectangle in Fig. 2(a)]. Fig. 12(a) is an enlarged view of the airport and vicinity. It is expected that buildings near an airport are low in height, including warehouses, repair factories, and some residential areas. A total of 8464 PS were within the rectangle. Three observations are noted. First, almost no PS were within the airport's runways and immediate vicinity, but PS representing the heights of surrounding buildings exist [see Fig. 12(b) and (c)]. Thus, the heights of the runway and the immediate areas cannot be derived. Second, although the heights derived by the PS-TSPA algorithm are mostly about zero in the area, as shown in Fig. 12(d), about 17% are negative. The negative heights are mainly concentrated at the edge of the airport. The negative height values jump to ∼89% after the Tomo-PS method. One possible reason is no PS are distributed in the airport, resulting in the arcs connecting points across the airport. Then, the arc lengths may be too long to estimate an accurate phase gradient. Moreover, a longer arc is more susceptible to the influence of decorrelation factors and atmospheric noise. Similar observations also exist on road surfaces and surrounding road areas, producing suboptimal or even unacceptable results in height estimation. Wu et al. [15] reported similar observations. Third, although the PS-TSPA algorithm may show a more robust tolerance for the long or "bad" arcs than the Tomo-PS, the ∼17% of the negative heights should still be too many. Thus, further studies are needed to reduce the negative values.

VI. CONCLUSION AND REMARKS
Combing the PS technique and the TSPA, we studied an algorithm for estimating building heights in urban areas with the time-series SAR observations of different baselines. The algorithm is called the PS-TSPA algorithm. It takes advantage of the global phase gradient information but does not consider the phase continuity assumption. Thus, the algorithm should be robust to the noise effect.
With the time series of 17 TerraSAR-X data, we estimated building heights in the southern urban area of Chengdu, China, where a wide range of urban targets (e.g., the railway station, road surfaces, overpasses, and low and high buildings) with different heights exist. The derived building heights were encouraging and satisfactory overall. In comparison, the Tomo-PS algorithm was used to estimate the heights in the same area and with the same SAR data. The study shows that although both methods estimated heights, significant height differences existed qualitatively and quantitatively.
With in situ measurements of four tall buildings and a railway station in the study area, we evaluated both methods in height estimations. The buildings and their surrounding environments vary from a simple and isolated building with minimum surroundings to a cluster of several buildings interplaying with surrounding buildings. In three out of four cases, the PS-TSPA algorithm outperformed the Tomo-PS method in the building height estimation. At the railway station, the PS-TSPA algorithm not only estimated the heights of the facility appropriately but also observed the height difference between the station roof and railway tracks, whereas the estimated heights after the Tomo-PS method were concerned. Thus, the former should be preferred in height estimation in urban areas.
Before ending, we would offer the following remarks. First, the backscatter from mono-state and side-looking SAR is none or weak when the surface is smooth or slightly rough because the incoming radar signal energy is mostly specularly reflected. Without adequate backscatter from the rooftop or ground surface, PS represent the height lacks. Then, the height estimation becomes unavailable. Second, the lack of PS within a large flat area can pose another challenge in height estimation. Without the scatterers, arcs of the network formed by the PS may be across the area. A long arc can negatively affect the phase gradient derivation. A long arc is also more susceptible to the decorrelation effect (e.g., the temporal variations and atmospheric conditions) in the time-series SAR data when unwrapping the phase [23]. To resolve the long arc effect, one needs to improve the Delaunay network algorithm further or evaluate and remove low-quality arcs [15], [23]. At the same time, the total unit modularity of the constraint matrix of the PU of PS should be fully achieved [40]. Finally, due to the computational limitation inhouse, solving the 16-absolute phase of PS could not be achieved simultaneously. Instead, we derived the scatterers using three groups, 6, 6, and 4 interferograms. The results were combined. The division might influence the algorithm's performance.