Transformative Urban Changes of Beijing in the Decade of the 2000s

The rapid economic growth, the exodus from rural to urban areas, and the associated extreme urban development that occurred in China in the decade of the 2000s have severely impacted the environment in Beijing, its vicinity, and beyond. This article presents an innovative approach for assessing mega-urban changes and their impact on the environment based on the use of decadal QuikSCAT (QSCAT) satellite data, acquired globally by the SeaWinds scatterometer over that period. The Dense Sampling Method (DSM) is applied to QSCAT data to obtain reliable annual infrastructure-based urban observations at a posting of ~1 km. The DSM-QSCAT data, along with different DSM-based change indices, were used to delineate the extent of the Beijing infrastructure-based urban area in each year between 2000 and 2009, and assess its development over time, enabling a physical quantification of its urbanization which reflects the implementation of various development policies during the same time period. Eventually, as a proxy for the impact of Beijing urbanization on the environment, the decadal trend of its infrastructure-based urbanization is compared with that of the corresponding tropospheric nitrogen dioxide (NO2) column densities as observed from the Global Ozone Monitoring Experiment (GOME) instrument aboard the second European Remote Sensing satellite (ERS-2) between 2000 and 2002, and from the SCanning Imaging Absorption SpectroMeter for Atmospheric CHartographY aboard of the ESA’s ENVIronmental SATellite (SCIAMACHY /ENVISAT) between 2003 and 2009. Results reveal a threefold increase of the yearly tropospheric NO2 column density within the Beijing infrastructure-based urban area extent in 2009, which had quadrupled since 2000.


Introduction
The Intergovernmental Panel on Climate Change (IPCC) Working Groups (WG) II and III recognize that a large fraction of greenhouse gases from power generation, industry, transportation, and consumption can be attributed to human settlements [1]. Due to the rapid growth of urban populations, and associated infrastructure, cities are becoming increasingly vulnerable to both climate change, which may exacerbate both natural and man-made hazards [2][3][4], and pollution [5,6].
The economic boom that took place in rapidly developing countries, such as those observed in China and India in the 2000s, and the associated great exodus from rural to urban areas, has led to a rapid and extreme urban development that has affected the environment not only in the urban vicinity, but also at regional [7,8] and, perhaps, global scales [9]. Nevertheless, major gaps exist in the assessment of the impact of urbanization on the environment at multiple spatial and temporal scales; mostly because of the difficulty of consistently estimating the extent and typology of urban areas through time and space [10] and the subsequent, and often biased, distribution of environmental monitoring stations [11,12].
Indeed, administrative extents of urban areas are often defined differently and inconsistently both within and among countries, which may not be representative of the underlying physical infrastructure [13]. Furthermore, the use of moderate spatial resolution optical remote sensing data, such as multi-temporal and globally available Landsat Enhanced Thematic Mapper Plus (ETM+) and Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER) images, to consistently estimate the type and extent of urban areas over space and time has proven to be a substantial challenge due to the presence of clouds, a recurrent phenomenon in tropical regions [14], the considerable intraurban and interurban spectral heterogeneity [15], and the spectral similarity among many impervious surfaces and their surrounding natural environments [16,17]. While high spatial resolution optical data can resolve individual components of the urban mosaic [18], such data are only available for small areas at irregular points in time [19]. Similarly, the National Oceanic and Atmospheric Administration (NOAA) Defense Meteorological Satellite Program Operational Line Scanner (DMSP-OLS) and the Suomi National Polar-Orbiting Partnership (NPP) Visible Infrared Imaging Radiometer Suite (VIIRS) night light data present a number of well-known limitations precluding their use to solely distinguish between different urban typologies and quantitatively assessing the horizontal and vertical growth of urban areas over space and time [20][21][22]. Synthetic aperture radar (SAR) data, which could potentially be used to consistently map infrastructure-based urban areas [23], are generally not collected globally over regular time intervals-with partial exceptions represented by the Shuttle Radar Topography Mission (SRTM) [24], which collected global SAR data in February 2000, and the TanDEM-X satellite, which collected global SAR data in 2011-2012 [25]. Very high-resolution spotlight SAR [26] and Lidar data [27] have been successfully used to map infrastructure-based urban areas at a meter or even sub-meter resolution, but their coverage is extremely limited in both space and time [27].
In this study, an innovative approach for consistently and quantitatively assessing infrastructure-based urbanization over space (i.e., horizontal and vertical urbanization/urban growth) and time, and its impact on the environment, is demonstrated for Beijing as a case study ( Figure 1). Such an approach is based on the use of (i) decadal level L1B Ku-Band (at the frequency of 13.4 GHz and wavelength of 2.24 cm) radar backscatter data, collected by the SeaWinds scatterometer aboard the QuikSCAT satellite (henceforth referred to as QSCAT data) [28], and (ii) decadal tropospheric nitrogen dioxide (NO 2 ) column densities, retrieved from the Global Ozone Monitoring Experiment (GOME) spectrometer on the European Remote Sensing Satellite 2 (ERS-2) [29,30] and the SCanning Imaging Absorption SpectroMeter for Atmospheric CHartographY (SCIAMACHY) aboard the European Space Agency's (ESA) ENVIronmental SATellite (ENVISAT) [31].
First, the Dense Sampling Method (DSM), developed at the Jet Propulsion Laboratory of the National Aeronautics and Space Administration [19], is applied to QSCAT slice data to obtain reliable annual physical infrastructure-based urban observations at 1-km posting. The DSM-processed QSCAT data (henceforth referred to as DSM data) are then used to (i) estimate the extent of the Beijing infrastructure-based urban area in each year of the 2000s, (ii) distinguish among different urban typologies, and (iii) define four different DSM-based indices enabling a quantitative assessment of land cover and intra-urban changes over time. Ultimately, to evaluate the impact of the Beijing infrastructure-based urbanization on the environment in terms of air pollution, its decadal urbanization trend is compared with that of the corresponding tropospheric nitrogen dioxide (NO 2 ) column densities. urbanization trend is compared with that of the corresponding tropospheric nitrogen dioxide (NO2) column densities.

Study Area
Beijing ( Figure 2) has experienced rapid economic growth in the past two decades (the city's gross domestic product increasing from 28.49 billion CNY in 1986 to 787 billion CNY at the end of 2006 [32]) and represents a prominent example of rapid urbanization in a strongly emerging economy. During the same period, Beijing's economic growth was accompanied by the development of its urban area, which was planned through three urban Master Plans in 1982, 1992, and 2004 [33].

Study Area
Beijing ( Figure 2) has experienced rapid economic growth in the past two decades (the city's gross domestic product increasing from 28.49 billion CNY in 1986 to 787 billion CNY at the end of 2006 [32]) and represents a prominent example of rapid urbanization in a strongly emerging economy. During the same period, Beijing's economic growth was accompanied by the development of its urban area, which was planned through three urban Master Plans in 1982, 1992, and 2004 [33]. urbanization trend is compared with that of the corresponding tropospheric nitrogen dioxide (NO2) column densities.

Study Area
Beijing ( Figure 2) has experienced rapid economic growth in the past two decades (the city's gross domestic product increasing from 28.49 billion CNY in 1986 to 787 billion CNY at the end of 2006 [32]) and represents a prominent example of rapid urbanization in a strongly emerging economy. During the same period, Beijing's economic growth was accompanied by the development of its urban area, which was planned through three urban Master Plans in 1982, 1992, and 2004 [33].  The city's development has been characterized by two primary phases that can be distinguished based on their urban growth rate and mechanism. In particular, the period 1996-2001 was characterized Remote Sens. 2020, 12, 652 4 of 24 by strong development along roads, followed by new city developments, and limited growth in areas around the city, while the period 2001-2006 experienced the highest development in riverside areas, followed by development in the city, and an increase of roadside built-up areas [35].
Furthermore, the urban development within the Beijing Municipality Area was significantly different from the that detailed in the plan. Indeed, during the two planning periods of 1983-1993 and 1993-2005, more urban development was carried out outside the planned urban growth boundaries than within [10]. Nowadays, Beijing Municipality covers a total area of 16,410 km2 with a population of about 21.5 million [36], which is almost twice that of 1976 [33]-this coincided with the number of vehicles increasing from 0.82 to more than 3.5 million since 1994 [37].
The simultaneous increase of vehicles and industrial activities has significantly increased the sources of emissions of both volatile organic compounds and NOx (NO+NO 2 ). Since 1998, Beijing adopted a series of measures in an effort to both control and decrease air pollution [32]. However, even though some air quality improvements were achieved by the end of the 1990s, air pollution worsened again during the 2000s due to the impressive growth rate that characterized the same period [7].
Beijing is also surrounded by three other highly populated, urbanized, and industrialized provinces, with four large cities having between 3.3 and 9.1 million inhabitants [7], which are considered significant contributors in the deterioration of Beijing's air quality [8,37,38].

QSCAT Data and Dense Sampling Method for Urban Observations
Level L1B Ku-Band QSCAT data, available through the NASA Jet Propulsion Laboratory Physical Oceanography Distributed Active Archive Center (PO.DAAC; https://podaac.jpl.nasa.gov/dataset/ QSCAT_LEVEL_1B_V2), were continuously collected between July 19th, 1999 and November 22nd, 2009 by the SeaWinds scatterometer.
The data were collected using a scanning pencil-beam antenna across a large swath (as wide as 1800 km) over Earth's surface, as the QuickSCAT satellite moved along drifting orbits around the Earth [39]. QSCAT data present a number of clear advantages, with respect to other remote sensing-based datasets, for providing consistent physical infrastructure-based urban observations without gaps in time or space [19,21]. Indeed, (i) QSCAT data have a daily coverage of about 90% of Earth's surface regardless of cloud coverage or level of light, (ii) QSCAT backscatter responses are dependent on the number, size, and construction materials of buildings and other structures (e.g., higher backscatter for more structures, larger and taller buildings, and stronger construction materials such as steel versus wood), and (iii) QSCAT drifting orbits allow the azimuth diversity in the backscatter measurements at each targeted area to detect the presence of buildings and other structures regardless of their directional alignment. Furthermore, QSCAT backscatter measurement is accurate to 0.2 dB (3-σ); [40], which is equivalent to approximately 1.57% in root-mean-square error, enabling QSCAT to detect large and rapid, as well as small and slow, environmental changes [41].
However, an important limitation of the raw QSCAT data is represented by their low spatial resolution-about 25 km in azimuth by 37 km in range [28], which is clearly inadequate to delineate and monitor inter-annual changes in intraurban areas. To overcome this limitation, DSM is used to increase the spatial resolution of the QSCAT data, at the expense of their temporal resolution (by reducing the latter from daily to annual), and thus to enable their use for assessing urbanization and its effect on multiple environmental matrices [13,[41][42][43][44].

of 24
Here, DSM is summarized to provide a basic understanding of the DSM data used in this study. At each location, with original QSCAT data densely compiled annually and posted at 30 arc seconds (approximately 1 km at the equator) in latitude and longitude (World Geodetic System 84; WGS84), a mathematical transform, called the Rosette Transform is applied to calculate an ensemble average of measured QSCAT backscatter values, (dB), using Equation (1): where N is the total number of scatterometer footprints centered at the considered location, and σ 0 (φ i , t i ) is the backscatter value measured at time t i and azimuth direction φ i of each footprint. In (1), σ 0M is the mean part obtained by the Rosette Transform, and represents the residual from the zero-mean fluctuating part of the radar backscatter. For a target like an urban area, having certain persistent characteristics over time whilst transient perturbations occur (e.g., new/less structures, traffic, rain, snowfall, etc.), σ 0M tends to a significant and stable value while becomes small with sufficiently large number of samples [19]. Thus, for each considered location, by sacrificing the temporal resolution of the QSCAT data over a long time period (over a year, for instance) in order to collect a sufficient number of samples, it is possible to (i) obtain a stable value for the mean part, and (ii) characterize the perturbations represented by the fluctuation term by an Index of Variability, IV, using Equation (2): where σ 0 (dB) is the linearized ensemble average of measured QSCAT backscatter values, and δ (dB) is the linearized standard deviation of σ 0 (dB) [19]. Thus, each DSM dataset is produced from a rectangular array of points, uniformly gridded at 30 arc seconds in latitude and in longitude (GCS-WGS84), with each point containing a stable σ 0 value and the associated IV. Here, QSCAT data acquired in the 2000s are used to generate annual DSM σ 0 and IV gridded datasets of the study area (Figures 3 and 4, respectively) bounded by latitudes 39.375 • S-40.5 • N and longitudes 115.75 • E-117.25 • E. In Figure 3 (and in subsequent figures), the "filtered" Beijing 1994-1995 stable night light-based extent [20], obtained using a 10% detection frequency threshold to include only persistent lights from the city while excluding gas flares and ephemeral events such as fires, is used as a common and consistent spatial reference. This is to better visualize the spatiotemporal change of the Beijing DSM-based urban area extent (see also Section 4.1). Refer to Figure S1, available in the Supplementary Materials, to visually compare annual DSM-based urban area extents against the corresponding intercalibrated stable night light data [45] for each year in the 2000s.
Beyond the capability of delineating boundaries demarcating DSM-based urban and rural/natural areas, σ 0 values can also be used to (i) distinguish among different intra-urban typologies, such as city cores (including commercial and industrial centers) with larger building volume and residential areas with lower building volume, and (ii) monitor their changes over time (see Section 4.2). This is because the value of σ 0 in urban areas is dependent on the 3D building volume, with a consistent linear relationship over a large dynamic range without a saturation problem [27]. It is noted that the capability of DSM data to capture physical urban patterns has been validated both in two-dimensional (2D) extent [46] and in three-dimensional (3D) building volume [27]. Remote Sens. 2019, 11, x FOR PEER REVIEW 6 of 25    [20] (black contour) is overlaid on all datasets as a common and consistent spatial reference area.

Tropospheric NO2 Columns
Tropospheric NO2 column densities can be measured from space using the Differential Optical Absorption Spectroscopy (DOAS) Method [47]. The data used in this study, representing tropospheric NO2 column densities averaged annually and posted at 7.5 arc minutes (approximately 14 km at the equator), are derived from GOME [29] for 2000-2002 and SCIAMACHY [31] for 2003-2009 [48]. GOME and SCIAMACHY tropospheric NO2 observations (expressed as 10 15 molec/cm 2 ) are available through the University of Bremen IUP DOAS data archive (http://doasbremen.de/no2_from_gome.htm and http://doas-bremen.de/no2_tropos_from_scia.htm, respectively) and can be successfully used to observe the temporal evolution of tropospheric NO2 columns worldwide [49][50][51]. It is important to note that although the uncertainty of absolute values for individual NO2 observations may be considerable, uncertainty is significantly reduced when averaging over a long period such as one year [52]. Furthermore, since many sources of uncertainty are systematic, relative values such as those used in this study, representing the variation of  [20] (black contour) is overlaid on all datasets as a common and consistent spatial reference area.

Tropospheric NO 2 Columns
Tropospheric NO 2 column densities can be measured from space using the Differential Optical Absorption Spectroscopy (DOAS) Method [47]. The data used in this study, representing tropospheric NO 2 column densities averaged annually and posted at 7.5 arc minutes (approximately 14 km at the equator), are derived from GOME [29] for 2000-2002 and SCIAMACHY [31] for 2003-2009 [48]. GOME and SCIAMACHY tropospheric NO 2 observations (expressed as 10 15 molec/cm 2 ) are available through the University of Bremen IUP DOAS data archive (http://doas-bremen.de/no2_from_ gome.htm and http://doas-bremen.de/no2_tropos_from_scia.htm, respectively) and can be successfully used to observe the temporal evolution of tropospheric NO 2 columns worldwide [49][50][51]. It is important to note that although the uncertainty of absolute values for individual NO 2 observations may be considerable, uncertainty is significantly reduced when averaging over a long period such as one year [52]. Furthermore, since many sources of uncertainty are systematic, relative values such as those Remote Sens. 2020, 12, 652 8 of 24 used in this study, representing the variation of tropospheric NO 2 column densities over time, have an even lower uncertainty, which is estimated to be around 15% over China [48].

Beijing Urban Extent
Although values of σ 0 are expected to be high within urban areas, where structures are located (e.g., houses, factories, shopping malls, skyscrapers, road network, etc.), in some cases it may be difficult to distinguish high-backscatter urban areas from high-backscatter natural areas such as forests or snowcapped mountain peaks [53]. Thus, backscatter values need to be combined with IV values to delineate the actual extent of infrastructure-based urban areas. Indeed, due to the anisotropy of structures, which behave as radar corner reflectors aligned in a certain azimuth direction, the angular variability of the backscatter is typically larger in urban areas than in natural environments where scatterer orientations are more random and scattering effects are mostly incoherent [19].
Hence, co-registered DSM σ 0 and IV gridded datasets can be combined to correctly identify grid cells characterized by the presence of structures. This can be done by identifying (i) an appropriate σ 0 threshold value to select all grid cells potentially representing urban areas (i.e., all grid cells characterized by backscatter values higher than the considered σ 0 threshold), and (ii) an IV threshold value to exclude cells with high-backscatter natural targets from previously selected all grid cells (i.e., to exclude all grid cells characterized by backscatter values higher than the considered σ 0 threshold but lower IV values than the considered IV threshold). In this study, to delineate the extent of the tropospheric NO2 column densities over time, have an even lower uncertainty, which is estimated to be around 15% over China [48].

Beijing Urban Extent
Although values of are expected to be high within urban areas, where structures are located (e.g., houses, factories, shopping malls, skyscrapers, road network, etc.), in some cases it may be difficult to distinguish high-backscatter urban areas from high-backscatter natural areas such as forests or snowcapped mountain peaks [53]. Thus, backscatter values need to be combined with IV values to delineate the actual extent of infrastructure-based urban areas. Indeed, due to the anisotropy of structures, which behave as radar corner reflectors aligned in a certain azimuth direction, the angular variability of the backscatter is typically larger in urban areas than in natural environments where scatterer orientations are more random and scattering effects are mostly incoherent [19].
Hence, co-registered DSM and IV gridded datasets can be combined to correctly identify grid cells characterized by the presence of structures. This can be done by identifying (i) an appropriate threshold value to select all grid cells potentially representing urban areas (i.e., all grid cells characterized by backscatter values higher than the considered threshold), and (ii) an IV threshold value to exclude cells with high-backscatter natural targets from previously selected all grid cells (i.e., to exclude all grid cells characterized by backscatter values higher than the considered threshold but lower IV values than the considered IV threshold). In this study, to delineate the extent of the Beijing urban area in each year from 2000 to 2009, each annual DSM gridded dataset (   For each year, the decadal (2000-2009) DSM IV gridded dataset, is used instead of the corresponding annual DSM IV dataset, to better distinguish between high backscatter grid cells representing urban areas and high backscatter grid cells representing natural/rural areas. By using stable decadal IV values, it is possible to identify a set of potential urban (S IV ) versus natural/rural grid cells for each Remote Sens. 2020, 12, 652 9 of 24 year, and then use the above-the-threshold σ 0 value to determine the actual urban grid cells in each annual S IV . It is noted that a natural/rural grid cell in 2000 not becoming urban will have a low decadal IV value, while a natural/rural grid cell becoming urban at some point will have a higher decadal IV value. On the contrary, the decadal IV value for a grid cell that was urban from the beginning will be high. Thus, the role of IV is to detect grid cells that are already urban or may become urban; for those pixels that can be urban, DSM σ 0 is used to monitor annual urban change throughout the 2000s.
To identify the optimal threshold values for σ 0 and decadal IV, both thresholds were iteratively adjusted by comparing the result from each iteration with high-resolution Landsat images of the city of Wuhan (China)-used as a test case for calibrating and validating the urban classification algorithm [27,43]. The annual results ( were compared to other independent estimates referring to these years. In particular, the 2000 DSM-based estimate is 6.3% bigger than the corresponding Landsat-based estimate of 1035 km 2 [55], the 2007 estimate is just 0.2% smaller than the Constructive Land area estimate of 3325.57 km 2 from the Beijing Municipal Bureau of Statistics [56], and the 2009 estimate is 7.6% smaller than the 2009 official urban area extent of 4480 km 2 that includes inner and  [20] and the Beijing administrative units (blue contour) [34] are overlaid on all datasets as a common and consistent spatial reference area.
The urban expansion was extensive in the city's southeastern sector and to a lesser degree in the northern districts of Huairou and Miyun. The presence of mountainous areas in the western and northwestern parts of the study area restricted the urban expansion in these directions.
The gridded DSM-based urban area extent datasets in Figure 6 and the ASPHAA algorithm [54] were then used to calculate the surface of the Beijing urban area extent in each year of the 2000s and its corresponding annual growth (Table 1). The surface of Beijing's DSM-based urban area extent quadrupled in the 2000s, with an increase of 3035 km 2 between 2000 and 2009, corresponding to an average annual increase of 337 km 2 . DSM-based estimates for 2000, 2007, and 2009 were compared to other independent estimates referring to these years. In particular, the 2000 DSM-based estimate is 6.3% bigger than the corresponding Landsat-based estimate of 1035 km 2 [55], the 2007 estimate is just 0.2% smaller than the Constructive Land area estimate of 3325.57 km 2 from the Beijing Municipal Bureau of Statistics [56], and the 2009 estimate is 7.6% smaller than the 2009 official urban area extent of 4480 km 2 that includes inner and outer suburbs-with the remaining 12,239.9 km 2 of the Beijing Municipality consisting of rural areas that include nature reserves and lightly populated mountains [57].
While the cumulative growth of the surface of the Beijing urban area extent has an overall linear trend in the 2000s (Figure 7), the year-to-year growth rate shows a large variability (Table 1)     In particular, the first major drop observed in 2005 (of about −12% with respect to that in the previous year) can be related to the growth reduction policy implemented in 2004 to mitigate the increasing social pressure and discontent associated with the mega rapid urbanization that was taking place in China [58]. The drop observed in 2007 (of about −4.5% with respect to that in the previous year) can be associated with the decision of the Chinese government to reassess the preparation of the 2008 Olympic Games, which led to the postponement of construction deadlines by a year [58]. Finally, the third major drop observed in 2009 (of about −14% with respect to that observed in the previous year), in turn, can be interpreted as the inherent consequence of the Olympic-related construction boom that occurred earlier in 2008. In particular, the first major drop observed in 2005 (of about −12% with respect to that in the previous year) can be related to the growth reduction policy implemented in 2004 to mitigate the increasing social pressure and discontent associated with the mega rapid urbanization that was taking place in China [58]. The drop observed in 2007 (of about −4.5% with respect to that in the previous year) can be associated with the decision of the Chinese government to reassess the preparation of the 2008 Olympic Games, which led to the postponement of construction deadlines by a year [58]. Finally, the third major drop observed in 2009 (of about −14% with respect to that observed in the previous year), in turn, can be interpreted as the inherent consequence of the Olympic-related construction boom that occurred earlier in 2008.

DSM-Based Change Indices
As mentioned in Section 3.1 and seen in Figure 3, for each year, the corresponding DSM σ 0 gridded dataset can be used to identify the inner core of Beijing (light-gray to red/light-orange), representing a well-established urban area with dense housing and tall buildings, its urban/suburban area (yellow to green) transitioning into its greater outskirt area with sparse houses (cyan), and the surrounding natural environment (blue). One can clearly observe the significant expansion and build-up characterizing the Beijing core with high backscatter [21] and the distinct concentric sprawl representing Beijing's urban development during the 2000s [10,59].
Thus, annual DSM σ 0 gridded datasets can also be used to monitor changes over time, both within urban areas and in their rural/natural surroundings. To this end, to quantitatively investigate these changes, four different DSM-based indices, described in the following subsections, were developed.  The Surface Change Index ∆σ 0 describes the total change which occurred between 2000 and 2009 (Figure 9). High positive values primarily indicate highly urbanized areas while negative values are likely associated with significant urban decay (e.g., demolition and/or abandonment) and/or environmental degradation (e.g., deforestation leading to denuded land). The ∆σ 0 gridded dataset (Figure 9)  The Surface Change Index ∆ describes the total change which occurred between 2000 and 2009 (Figure 9). High positive values primarily indicate highly urbanized areas while negative values are likely associated with significant urban decay (e.g., demolition and/or abandonment) and/or environmental degradation (e.g., deforestation leading to denuded land). The ∆ gridded dataset (Figure 9) [20]. The + sign marks the "hotspot" location showed in Figure 10.  [20]. The + sign marks the "hotspot" location showed in Figure 10.
As an important capability, the ∆ gridded dataset can be used to identify urbanization "hotspots", characterized by high values of ∆ , where the largest changes between 2000 and 2009 occurred due to significant build-up. One of these urbanization "hotspots", where many large buildings were constructed between 2003 and 2009, is detected in the westernmost part of the province of Tongzhou (plus sign in Figure 9) and verified with high-resolution images ( Figure 10). As an important capability, the ∆σ 0 gridded dataset can be used to identify urbanization "hotspots", characterized by high values of ∆σ 0 , where the largest changes between 2000 and 2009 occurred due to significant build-up. One of these urbanization "hotspots", where many large buildings were constructed between 2003 and 2009, is detected in the westernmost part of the province of Tongzhou (plus sign in Figure 9) and verified with high-resolution images ( Figure 10).
The second index, computed at the grid cell level and expressed in decibel per year (dB/year), is calculated as the slope of the linear regression line obtained by plotting annual DSM ∆σ 0 grid cell values against time over all years in the 2000s. This index is calculated using Equation (4): where d[F(σ 0 , t)]/dt represents the time derivative of the linear regression function F(σ 0 , t) of σ¯0 in each grid cell. The Surface Change Index S(σ 0 ) describes the rate of change of the land cover throughout the 2000s. Large positive slope values primarily indicate urban development within the corresponding grid cells, with greater slope values associated with significantly higher urban development. In contrast, null or negative slope values indicate grid cells typically associated with the steadier rural and natural areas surrounding Beijing.
Since S(σ 0 ) employs all annual data in the 2000s, rather than just the difference between two years as in (3), it has the advantage of being more stable compared to the noisier results of ∆σ 0 . This is evident when looking at the urban "hotspot" located in the Tongzhou province of Beijing, where the S(σ 0 ) gridded dataset contains a more consistent set of very high value grid cells (i.e., red cells in Figure 11) compared to that observed in the ∆σ 0 gridded dataset, which shows an inconsistent number of very high value grid cells (i.e., red cells) mixed with others having lower values (Figure 9 vs. Figure 11). This result illustrates the advantage of DSM data providing observations without gaps in time and space for rural/urban monitoring. The coefficient of determination R 2 associated with the regression line is calculated as an additional data product from the ( ) regression. Since R 2 represents how well the response variable variation can be explained by a linear model, high values of R 2 indicate a consistent linear trend that can be associated with planned urban development. Conversely, low values of R 2 indicate more random natural characteristics, with transient or intermittent variations, in areas remaining undeveloped. The R 2 gridded dataset ( Figure 12) seems to confirm that the urban development in The coefficient of determination R 2 associated with the regression line is calculated as an additional data product from the S(σ 0 ) regression. Since R 2 represents how well the response variable variation can be explained by a linear model, high values of R 2 indicate a consistent linear trend that can be associated with planned urban development. Conversely, low values of R 2 indicate more random natural characteristics, with transient or intermittent variations, in areas remaining undeveloped. The R 2 gridded dataset ( Figure 12) seems to confirm that the urban development in both the inner core of Beijing and its surroundings are the results of structural build-up.  [19], UDI(y) represents the total urban growth/development integrated over the existing urban area extent, with more and larger buildings, as well as over the nonurban areas that would become urban by 2009. The persistently increasing UDI(y) for Beijing from 2000 to 2009, that closely follows a linear trend (Figure 13), indicates that in this period Beijing experienced a continuous constant urban growth, in terms of both number and size of its physical structures. The drop observed in 2005 ( Figure  8) represents an exceptional case, as a consequence of a temporary policy change in the overall decadal urban build-up and expansion of Beijing [58].  [19], UDI(y) represents the total urban growth/development integrated over the existing urban area extent, with more and larger buildings, as well as over the nonurban areas that would become urban by 2009. The persistently increasing UDI(y) for Beijing from 2000 to 2009, that closely follows a linear trend (Figure 13), indicates that in this period Beijing experienced a continuous constant urban growth, in terms of both number and size of its physical structures. The drop observed in 2005 (Figure 8) represents an exceptional case, as a consequence of a temporary policy change in the overall decadal urban build-up and expansion of Beijing [58]. Second, the Building Development Index (BDI), for each year in the 2000s, represents a quantitative indicator of having more and/or larger structures, being built within the DSM-based urban area extent delineated in each year (see Section 4.1). The yearly BDI is calculated by averaging, for a given year, the values of all DSM grid cells located within the urban area extent referring to the same year. This index is calculated using Equation (6): where ( ) is the DSM backscatter value (dB) of a grid cell i classified as urban in the year y and N(y) is the total number of grid cells classified as urban in the year y. As defined, BDI(y) accounts only for what occurs within the urban area extent occupied by buildings and other structures in each year, while excluding areas classified as nonurban. Thus, as expected, in each year the value of BDI (Figure 14) is higher than the corresponding UDI value ( Figure  13), since the latter accounts for both urban and nonurban areas. Therefore, BDI(y) tracks the vertical built-up areas more closely than UDI(y), with the latter being more representative of the general longterm urban development. This is evident in that the BDI values in 2005 and 2006 (see ellipse in Figure  14) are similar as a direct consequence of the curb on urban development in 2005 [58].
The results from the combination of UDI and BDI highlight two simultaneous processes characterizing Beijing's urban development in the 2000s including (i) the urban expansion/sprawl due to the transformation of nonurban to urban areas and (ii) the densification/upward height of buildings due to increasing construction activities in already existing urban areas. Second, the Building Development Index (BDI), for each year in the 2000s, represents a quantitative indicator of having more and/or larger structures, being built within the DSM-based urban area extent delineated in each year (see Section 4.1). The yearly BDI is calculated by averaging, for a given year, the values of all DSM σ 0 grid cells located within the urban area extent referring to the same year. This index is calculated using Equation (6): where σ 0 i (y) is the DSM backscatter value (dB) of a grid cell i classified as urban in the year y and N(y) is the total number of grid cells classified as urban in the year y. As defined, BDI(y) accounts only for what occurs within the urban area extent occupied by buildings and other structures in each year, while excluding areas classified as nonurban. Thus, as expected, in each year the value of BDI (Figure 14) is higher than the corresponding UDI value (Figure 13), since the latter accounts for both urban and nonurban areas. Therefore, BDI(y) tracks the vertical built-up areas more closely than UDI(y), with the latter being more representative of the general long-term urban development. This is evident in that the BDI values in 2005 and 2006 (see ellipse in Figure 14) are similar as a direct consequence of the curb on urban development in 2005 [58]. Confined within the urban extent of 2000, here referred to as the inner core of Beijing, VDI characterizes its physical growth, between 2000 and 2009, due to more and larger structures. As such, VDI(y) specifically tracks the mostly vertical build-up, represented by an increase in the total building volume, within the inner core of Beijing ( Figure 15).
Of all of the development indices described above, VDI(y) has that largest value in each year of the 2000s, confirming that the internal core of Beijing contains the most and the largest build-up. Moreover, VDI(y) has the largest rate of change (i.e., slope of the regression fit function), compared to those of UDI and BDI, revealing the continuous and strong development of the inner core of Beijing. Together, all of the development indices contribute to a more complete characterization of Beijing urban change in the 2000s, from the Beijing inner core to its outward expansion. The results from the combination of UDI and BDI highlight two simultaneous processes characterizing Beijing's urban development in the 2000s including (i) the urban expansion/sprawl due to the transformation of nonurban to urban areas and (ii) the densification/upward height of buildings due to increasing construction activities in already existing urban areas.
Third, the Vertical Development Index (VDI), calculated for each year in the 2000s with respect to year 2000, represents a quantitative urbanization measure accounting for the vertical growth within the inner core of Beijing identified as urban in 2000. The yearly VDI(y) is calculated by averaging all DSM σ 0 grid cell values within the 2000 Beijing urban area extent using Equation (7): where σ 0 i (y) is the DSM backscatter value (dB) of a grid cell i in the year y classified as urban in 2000 and N(2000) is the total number of grid cells classified as urban in 2000. Confined within the urban extent of 2000, here referred to as the inner core of Beijing, VDI characterizes its physical growth, between 2000 and 2009, due to more and larger structures. As such, VDI(y) specifically tracks the mostly vertical build-up, represented by an increase in the total building volume, within the inner core of Beijing ( Figure 15).
Of all of the development indices described above, VDI(y) has that largest value in each year of the 2000s, confirming that the internal core of Beijing contains the most and the largest build-up. Moreover, VDI(y) has the largest rate of change (i.e., slope of the regression fit function), compared to those of UDI and BDI, revealing the continuous and strong development of the inner core of Beijing. Together, all of the development indices contribute to a more complete characterization of Beijing urban change in the 2000s, from the Beijing inner core to its outward expansion.

Beijing NO2 Pollution
Rapid urbanization and economic growth, along with the associated NOx emissions, have a key role in the increase of tropospheric NO2 [47,49,50]. High concentrations of NO2 in the troposphere are harmful to human health by both aggravating respiratory diseases in the short term and contributing to onset of asthma in the long term [60].
For each year in the 2000s, Figure 16 presents satellite observations of yearly tropospheric NO2 column densities (see Section 3.2) over the study area, overlapped by the corresponding Beijing DMSbased urban area extent (see Section4.1). Tropospheric

Beijing NO 2 Pollution
Rapid urbanization and economic growth, along with the associated NOx emissions, have a key role in the increase of tropospheric NO 2 [47,49,50]. High concentrations of NO 2 in the troposphere are harmful to human health by both aggravating respiratory diseases in the short term and contributing to onset of asthma in the long term [60].
For each year in the 2000s, Figure 16 presents satellite observations of yearly tropospheric NO 2 column densities (see Section 3.2) over the study area, overlapped by the corresponding Beijing DMS-based urban area extent (see Section 4.1).
Tropospheric NO 2 column density values were low in the early 2000s while increasing markedly between 2003 and 2009 ( Figure 16), along with Beijing's urban development during the same period. Nevertheless, an obvious decrease can be observed in 2008 (Figure 16), corresponding to the year of the Olympic Games. Indeed, because of the high levels of air pollution in Beijing prior to 2008, a series of control measures were implemented to substantially reduce air pollution during the pre-Olympic and the Olympic period [61,62]. In 2009, however, right after the Olympic year, tropospheric NO 2 column density values starkly increased again to a more serious level compared to 2007 ( Figure 16).
As a quantitative assessment of tropospheric NO 2 pollution, column densities (10 15  A regression analysis between yearly UDI values and yearly tropospheric NO2 columns, both calculated with respect to the Beijing DSM-based urban extent in 2009, is carried out to investigate the relation between Beijing urban development and tropospheric NO2 pollution ( Figure 18).  Over the 2000s, regression results indicate an increasing trend in tropospheric NO2 pollution alongside the urbanization of Beijing, with a high correlation coefficient of 0.9063. The R 2 value of 0.821 explains 82.1% of the tropospheric NO2 column trend with respect to the Beijing urban development trend.

Conclusions
Mega-urban changes in Beijing and their impacts on the environment are assessed consistently and quantitatively through an innovative approach using decadal QuickSCAT data and GOME/SCIAMACHY-derived tropospheric nitrogen dioxide (NO2) columns. Annual observations  Over the 2000s, regression results indicate an increasing trend in tropospheric NO2 pollution alongside the urbanization of Beijing, with a high correlation coefficient of 0.9063. The R 2 value of 0.821 explains 82.1% of the tropospheric NO2 column trend with respect to the Beijing urban development trend.

Conclusions
Mega-urban changes in Beijing and their impacts on the environment are assessed consistently and quantitatively through an innovative approach using decadal QuickSCAT data and GOME/SCIAMACHY-derived tropospheric nitrogen dioxide (NO2) columns. Annual observations

Conclusions
Mega-urban changes in Beijing and their impacts on the environment are assessed consistently and quantitatively through an innovative approach using decadal QuickSCAT data and GOME/SCIAMACHY-derived tropospheric nitrogen dioxide (NO 2 ) columns. Annual observations of Beijing urban area extent are based on the Dense Sampling Method (DSM) that allows a major increase in the QSCAT data sampling, from about 25 km 2 to a posting of about 1 km 2 , which is appropriate to map urban extent based on the presence of buildings and other structures rather than on the administratively defined city boundaries.
Compared to independent sources in multiple years, including results obtained using a Landsat-based method and government statistics, DSM applied to QuickSCAT data from the 2000s enables us to: (i) Monitor the annual growth of Beijing's urban area extent, which quadrupled in the 2000s with an increase of 3035 km 2 between 2000 and 2009, corresponding to an average annual increase of 337 km 2 . DSM-based results reflect the Chinese political and policy changes affecting Beijing's urban development over the decade. These include the growth reduction policy implemented in 2005 and the pollution-curbing policy adopted before the Games of the XXIX Olympiad in 2008.
(ii) Observe the significant urban expansion and built-up in Beijing, including its distinct concentric sprawl during the 2000s. DSM backscatter values, combined with the associated Index of Variability values, enable us to distinguish between the inner core of Beijing (representing a well-established urban area with dense housing and tall buildings), the urban and suburban areas transitioning into Beijing's greater outskirt area (characterized by the presence of sparse houses), and the surrounding natural environment.
(iii) Infer processes of urban development (e.g., expansion, densification, and vertical growth) spanning both the inner core of Beijing and its surrounding areas using DSM-based development indices (i.e., Urban Development Index, Building Development Index, and Vertical Development Index). Results show large changes both in the inner urban core and urbanized areas of Beijing, while suggesting that the rural/natural areas around Beijing remained relatively stable in the 2000-2009 time period. Given that DSM can be used to process the whole global QuikSCAT data collected in the 2000s, similar analyses could provide more details and insights on the urbanization rate and type of urban development for the whole of China and other megacities around the world.
Furthermore, the DSM-based urban analysis allows a consistent comparison with the year-by-year changes of tropospheric NO 2 column densities measured from space. This study highlights a clear relationship between the Beijing urban growth/development and its tropospheric NO 2 pollution trend. The comparison reveals a triple increase in the yearly tropospheric NO 2 column densities spatially averaged over the 2009 Beijing DSM-based urban extent, which quadrupled since 2000. Furthermore, the results capture the effects of the pollution-curbing policy implemented before the Olympic year (i.e., 2008) as well as the sharp increase in tropospheric NO 2 pollution in 2009, with the latter clearly indicating that the series of measures developed to control and decrease Beijing air pollution in preparation for and during the Games of the XXIX Olympiad were temporary. This analysis on the relation between tropospheric NO 2 pollution and urban growth/development highlights the capability of multi-source imaging from space to quantitatively assess urbanization impacts on the environment. Future research will be necessary to address several limitations in current methods to (i) capture the urban-rural transition as a gradient rather that an abrupt binary boundary while using the Index of Variability to explicitly quantify the larger deviation from the mean value for more urbanized areas compared to that of more rural and natural areas, (ii) delineate fine features and their change in urban areas with the availability of high-resolution products such as the Global Urban Footprint [25] and time-series remote sensing data such as Sentinel [63], and (iii) extent long-term urban monitoring with future satellite missions such as the NASA-Indian Space Research Organisation (ISRO) Synthetic Aperture Radar (SAR) Mission [64]. Nevertheless, the use of DSM-processed QSCAT data offers a significant step toward a consistent and quantitative assessment of both horizontal and vertical urbanization that can be used to complement other currently available remote sensing-based urban and settlement datasets [21] obtained using different approaches and input data [25,[65][66][67]. Multiple current and future satellite X-band SAR missions, such as the TerraSAR-X (launched), TanDEM-X, COnstellation of small Satellites for Mediterranean basin Observation (COSMO-SkyMed), Advanced Satellite with New system Architecture for Observation-2 (ASNARO-2), and LOTUSat-1 and LOTUSat-2, will provide radar backscatter global data useful for 3D urban observations at a resolution of 10-100 meters [27].