Abstract
Aerodynamic roughness is very important to urban meteorological and climate studies. Radar remote sensing is considered to be an effective means for aerodynamic roughness retrieval because radar backscattering is sensitive to the surface roughness and geometric structure of a given target. In this paper, a methodology for aerodynamic roughness length estimation using SAR data in urban areas is introduced. The scale and orientation characteristics of backscattering of various targets in urban areas were firstly extracted and analyzed, which showed great potential of SAR data for urban roughness elements characterization. Then the ground truth aerodynamic roughness was calculated from wind gradient data acquired by the meteorological tower using fitting and iterative method. And then the optimal dimension of the upwind sector for the aerodynamic roughness calculation was determined through a correlation analysis between backscattering extracted from SAR data at various upwind sector areas and the aerodynamic roughness calculated from the meteorological tower data. Finally a quantitative relationship was set up to retrieve the aerodynamic roughness length from SAR data. Experiments based on ALOS PALSAR and COSMO-SkyMed data from 2006 to 2011 prove that the proposed methodology can provide accurate roughness length estimations for the spatial and temporal analysis of urban surface.
1. Introduction
Urban areas occupy only a very small fraction of the Earth’s land area, yet 60–80% of final energy and concentrate materials, wealth, and innovation is consumed in urban areas. With rapid urbanization, the population of urban areas has exceeded rural areas in the past decade [1]. The increased amount of roughness elements, which are mainly composed of buildings and road networks in urban areas, dramatically alters the surface roughness and urban boundary layer dynamics [2], which in turn results in more frequent air pollution and extreme weather. More specifically, the increased drag and enhanced turbulent momentum flux result in reduced wind speeds within the urban boundary layer. The drag effect of roughness elements is very important for urban mesoscale numerical models and it is usually parameterized using an aerodynamic roughness length. Aerodynamic roughness length is usually defined as the height where the wind velocity is equal to zero. It is an important aerodynamic parameter and reveals the exchange between the atmosphere and land surfaces. Precise mapping of aerodynamic roughness in urban areas is crucial for different modeling applications, like urban climate or air quality studies, as well as for planning purposes such as the identification of potential ventilation paths [3]. Nevertheless, urban areas are the roughest of all aerodynamic boundaries, and the spatial distribution of roughness elements in urban areas is highly heterogeneous [4]. It is very difficult to quantitatively retrieve the aerodynamic roughness in urban areas and its study is in a very preliminary stage [5].
Most of the conventional methods used in aerodynamic roughness calculations are meteorological methods, which utilize the wind and temperature gradient data acquired by flux towers or meteorological stations. The cost of this approach is great, and the calculated aerodynamic roughness length only represents a very small extent of an entire urban area. Other kinds of methods calculate aerodynamic roughness using morphological methods according to the height, distribution, and frontal areas of the roughness elements. Nevertheless, most current methods assume that the height of roughness elements is the same and they are distributed evenly, which is very different to the actual high degree of heterogeneity that is presented in urban areas. Infrastructure spatial data and LiDAR data are anticipated to overcome this limitation, yet this kind of data is expensive and its update period is very long, which restricts its application within very small areas. As such, there is an urgency to find an effective method for the aerodynamic roughness mapping in urban areas [6].
Synthetic Aperture Radar (SAR) is anticipated to provide a valuable tool for roughness mapping because the backscattering received by the SAR sensor is sensitive to the geometric structure of the target, which in turn is directly related to surface roughness. And SAR data has been used for aerodynamic roughness retrieval of uniform surfaces, such as deserts and vegetation areas [5, 7]. Yet aerodynamic roughness estimation in urban areas is very difficult because of the high heterogeneity of roughness elements, and the challenges related to high heterogeneity of urban areas still need to be resolved. The current methods for urban roughness mapping can be classified into three categories.
The first category maps the aerodynamic roughness through first classifying SAR images and then specifying roughness values for each class according to a look-up table. One of the most famous methods is Davenport method [8], yet the 8 classes can not precisely describe the roughness elements in urban areas. Zoran used multiple SAR images to map urban aerodynamic roughness through precisely identifying non-built-up areas, residential areas, industry areas, and large groups of buildings [9]. Stetson estimated the aerodynamic roughness length in Houston metropolitan area through applying a coefficient to the commonly used roughness length of a given class, which was determined according to the quotient of a pixel’s backscattering and the mean backscattering of the corresponding land use type [10]. Bidaut et al. used multitemporal SAR images for supervised classification and then mapped the urban aerodynamic roughness using a look-up table for each land use class [11]. This is the most widely used method, yet it is not suitable for mesoscale models because urban areas need more precise classification results and yet it is very difficult because urban areas are highly heterogeneous.
The second category calculates aerodynamic roughness using a morphological method that utilizes the geometric parameters of roughness elements such as the height, plan area index, and front area index. Basly et al. mapped the aerodynamic roughness of Nantes using ERS images through extracting structural features and land use class attributes [12]. Bechtel et al. used interferometric SAR data to extract the topology, statistics, and texture of roughness elements and then used the mean and the variance of the height of the roughness elements within a sector to map the aerodynamic roughness in urban areas [13], and Langkamp further analyzed the effect of various dimension of sector area on aerodynamic roughness retrieval [14]. This kind of method exerts strict requirements for the precision of interferometric measurements which is generally regarded very difficult for urban areas.
The last category tries to build a quantitative relationship between aerodynamic roughness and SAR backscattering. Jeyachandran et al. first analyzed the correlation between SAR backscattering and the urban canopy parameters such as the mean height, plan area fraction, and frontal area index of roughness elements in urban areas and then built exponential relationship model between the SAR backscattering and the aerodynamic roughness length [15]. And experiments showed that this method could provide roughness parameters acceptable for meteorological models. Jeyachandran et al. defined grid cells ranging in size from 25 to 250 m for correlation analysis between SAR data and urban roughness but did not consider aerodynamic roughness in heterogeneous urban areas is highly related to the roughness elements in the upwind area and the scale and orientation sensitivity analysis is indispensable.
Aimed at addressing these limitations, this paper proposed an aerodynamic roughness estimation methodology from SAR data based on a precise scale and orientation analysis, where ground truth aerodynamic roughness values were calculated using meteorological tower data. First, the scale and orientation characteristics of backscattering of roughness elements in urban areas were thoroughly analyzed. Then, the optimal dimension of the upwind sector for aerodynamic roughness elements calculation was determined through a correlation analysis of 22 scenes of ALOS PALSAR data obtained from 2006 to 2011 and the wind gradient data acquired by the Beijing meteorological tower. Finally a quantitative relationship between SAR backscattering and aerodynamic roughness was set up to map the spatial and temporal aerodynamic roughness length of Beijing city from ALOS PALSAR and COSMO-SkyMed data. This paper is organized as follows. In Section 2, we introduce the study area, the SAR data, and the wind gradient data acquired from the Beijing meteorological tower. Then in Section 3, the principles and the proposed methodology for aerodynamic roughness estimation from SAR data are introduced. Experimental results and a discussion are presented in Section 4, and our conclusions are given in Section 5.
2. Study Area and Experimental Data
2.1. Study Area
Beijing city is the capital of China and one of the main megacities in the world. Since the reform and open policy was implemented in China, land use and land cover of Beijing city have changed dramatically due to rapid social and economic development. Buildings have become increasingly taller, and the spatial distribution of infrastructures in Beijing city becomes more and more divergent. These changes have dramatically altered the surface roughness and thus induced an obvious change in the urban boundary layer dynamics, such as changes in wind and temperature, which contributes to the deterioration in the air quality and frequent meteorological disasters in Beijing. The 325 m high Beijing meteorological tower has provided a powerful tool for meteorological and climate studies in Beijing. In this paper, wind gradient data observed by the meteorological tower was used to calculate the ground truth aerodynamic roughness. Figure 1(a) shows the Gaofen image of Beijing city on 12 October 2015, and (b) shows the area around the Beijing meteorological tower (14 km × 14 km) which is used for the correlation and quantitative relation analysis.
2.2. Wind Gradient Data
In Figure 1 the white star illustrates the position of the Beijing meteorological tower, where its longitude is 116°22′E and latitude is 39°58′N. This tower has made great contributions to meteorological studies of Beijing city. At the early stage of its establishment, this meteorological tower was in the suburb of Beijing city, and most targets around this tower were crops or bare soil. With the development of Beijing city, the uniform surface around the tower has changed into a very complicated and heterogeneous one. Figure 2 shows photos of the underlying surface in the northeast (0–45°) and southwest (180–225°) directions, respectively, which shows that there are great differences between the roughness elements in different directions. Yet the aerodynamic roughness value calculated using wind gradient data can only represent the aerodynamic roughness characteristics of the limited extent around the meteorological tower, and only using one value can not reveal the detailed information of the whole Beijing city due to the high spatial variations of the urban targets.
(a) Photo in the orientation of 0–45°
(b) Photo in the orientation of 180–225°
Figure 3 shows photos of the Beijing meteorological tower, which includes 15 observation layers with height, respectively, at 9 m, 15 m, 32 m, 47 m, 65 m, 80 m, 103 m, 120 m, 140 m, 160 m, 180 m, 200 m, 240 m, 280 m, and 320 m. At each layer, the anemometer, anemoscope, thermometer, and hygrometer are mounted. The sampling time interval for meteorological data is 20 seconds, and the observed parameters include wind velocity and direction, humidity, and temperature. In this paper, multilayer wind data five days before and after the SAR imaging date (totally 106 days for ALOS PALSAR and 220 days for COSMO-SkyMed) were acquired and used for the ground truth aerodynamic roughness calculation.
(a) The Beijing meteorological tower
(b) The wind measurement sensors
2.3. SAR Data
In this paper, a total of 22 scenes of Single Look Complex (SLC) ALOS PALSAR data collected from 2006 to 2011 were used in our analysis and modeling. All of the ALOS PALSAR data are in the L-band, StripMap mode, and the incidence angle is 38.7°. Among them, 10 scenes are single polarization mode (HH), and the other 12 scenes are dual polarization mode (HH and HV). The detailed parameters of the ALOS PALSAR data are shown in Table 1.
Besides ALOS PALSAR data, a total of 22 scenes of COSMO-SkyMed data from 2009 through 2011 were jointly used for the temporal variation analysis of the aerodynamic roughness of Beijing city. All of the COSMO-SkyMed data are in the X band, HH polarization, StripMap mode, and the incidence angle is 18.3°. The detailed parameters of the COSMO-SkyMed data are shown in Table 2.
3. Principles and Methodology
3.1. Basic Principles
Radar backscattering received by SAR sensor is determined by the geometric structure, dielectric characteristics of targets, as well as the sensor parameters, and the spatial relations between the SAR sensor and targets. And there are great differences between the backscattering mechanisms of various land use types. For smooth surfaces, such as water or smooth roads, specular reflection is the dominant scattering. For rough surfaces, such as trees or crops, volume scattering is the dominant scattering. And for very rough surfaces, such as built-up areas or other artificial targets, dihedral scattering is the dominant scattering. Besides, the backscattering characteristics of urban areas highly rely on the spatial scale. On low or middle resolution SAR image, the backscattering coefficient reflects the general characteristics of land use types, while, on high resolution SAR image, more detailed information can be revealed. Jeyachandran et al. used SIR-C data and analyzed the backscattering coefficients of typical targets in urban areas in 2009, and their results showed that the backscattering coefficient increases with the roughness of targets [15].
We extracted backscattering coefficients of typical land use types in Beijing city from ALOS PALSAR data and COSMO-SkyMed data, including residential areas, central business district, high-rise building area, green land, road, and water. And the high resolution optical data and field surveys were used for validation. Table 3 illustrates the results extracted from ALOS PALSAR data. From Table 3 we can see the mean values of the six land use types are obviously different, but the minimum and maximum values of some land use types are approaching. In particular the maximum values of the first three classes, that is, residential area, central business district, and high-rise building area, are the same. This is because dihedral scattering, trihedral scattering, and multiple bounce scattering between adjacent urban targets make the backscattered echoes superimposed together and cause saturation in some cases. This is easy to happen for tall building areas, so the maximum backscattering values of the first 3 classes are the same, which means the backscattered echoes are saturated. Therefore we think the mean value of the backscattering coefficient can better represent the characteristics of land use types in urban areas and reveal their difference. The average backscattering coefficient of high-rise building areas is the largest, with a mean value of 0.07 dB. The backscattering coefficient of the central business district has the second largest value, with a mean of −0.53 dB, and the backscattering coefficient of an ordinary residential area has the third largest value, with a mean of −6.23 dB. There are obvious differences among the average backscattering coefficients of various built-up areas, and they are all much higher than the backscattering coefficient of relatively smooth surfaces such as roads and waters.
Experiments based on COSMO-SkyMed showed that the absolute value of backscattering coefficient was a little different to that from ALOS PALSAR, yet the rule of the backscattering coefficient difference between various land use types was the same. Our results also were compared with that of Jeyachandran et al., which showed good accordance between them. That is to say, the average backscattering coefficients of SAR data can reveal the difference of various land use types in built-up areas. Thus, we think SAR data has great potentials for being used in aerodynamic roughness length retrieval because both the backscattering and the aerodynamic roughness are closely related to the geometric structures of roughness elements presented in a given land surface.
3.2. Aerodynamic Roughness Length Estimation Methodology from ALOS PALSAR Data
The procedures for aerodynamic roughness length estimation from SAR data are illustrated in Figure 4, which mainly include ALOS PALSAR data processing and backscattering analysis, ground truth aerodynamic roughness length calculation using the wind gradient data, correlation analysis between the backscattering coefficient and the ground truth aerodynamic roughness length for different upwind sectors of various dimensions, and ultimately establishing of a quantitative relationship and the aerodynamic roughness length mapping. The details are described as follows.
Processing ALOS PALSAR data: the multilook processing was first applied to the SLC ALOS PALSAR data to produce intensity images with an improved radiometric resolution and reduced speckle noises through averaging along the range and azimuth direction. When calculating the number of looks for multilook processing, the range resolution is firstly calculated using pixel spacing in slant range/sin (incidence angle). For the single polarized data, taking the image acquired on 23 June 2006 as an example, the range resolution is equal to 4.68/sin m. The number of looks in range direction is set to be 1, and the azimuth resolution should be the same with the range direction to obtain a square pixel size, that is, . The number of looks can be calculated in the similar way and then the multilook operation can be performed for the dual polarization mode data. So the multilook factors in range and azimuth for the single polarization images are one and two, respectively, and the multilook factors in range and azimuth for the dual polarization images are one and five, respectively.
Then, the 22 scenes of ALOS PALSAR images were coregistered, which required spatial registration and resampling (in cases where the pixel sizes differed) to correct relative translational shift, rotation, and scale differences. Then the De Grandi filtering was performed to these 22 scenes of multitemporal ALOS PALSAR images to exploit the space-varying temporal correlation of speckles between the images to further suppress the system-inherent multiplicative noise presented in the SAR data [16]. Finally, a geometric correction and radiometric calibration were applied to the ALOS PALSAR images to make a meaningful intercomparison of the backscattering coefficient determined at different dates. In this paper we did not use the HV channel data, and the backscattering coefficient of HH polarization was used for all analysis.
Processing of the wind gradient data observed by the meteorological tower: first, obvious error data was eliminated, and then the Richardson Number was used to determine atmosphere stability. Then gradient wind data classified as arising from neutral conditions were adopted to calculate the ground truth aerodynamic roughness length at different orientations using fitting iteration method [17].
According to Monin-Obukhov similarity theory, a logarithmic relationship exists between horizontal wind speeds and measurement height; that is,where is the friction velocity determined from the wind profile measurement (m/s), is von Karman’s constant (usually set to 0.4), is the mean wind speed (m/s) at height (m), is the zero-plane displacement (m), and is the roughness length (m).
Equation (1) can be rewritten aswhere there is usually , and takes 5~7.
Thus and can be calculated through two steps according to (2) using a fitting iteration method [18]. Firstly, supposing various values for , a linear regression analysis was performed between and , and the value of was taken as optimal when the maximum correlation between and was achieved. Then once the value of was decided, were calculated according to the slope of the linear regression line described as , and were derived from the intercept of the linear regression line described as , respectively. In this paper, the observed wind data at meteorological tower heights of 32 m, 47 m, 65 m, and 80 m were used for calculating the ground truth roughness length because they could more easily satisfy the logarithmic relationship [17].
The aerodynamic roughness is related to the topology of the roughness elements in the upwind sector, which can be described by three parameters: the central angle , the radius , and the opening angle [13, 19]. The dimension of the upwind sector also represents the effective range of the upwind, meaning out of which the local flow effects can be neglected. Zhu et al. suggested that the dimension of the upwind sector should be determined according to the height of the mounted anemometer and suggested a radius that is 100 times greater than this height, with an opening angle of 30° [19]. Bechtel et al. suggested that the dimension of the upwind sector should be determined according to the height of roughness elements and suggested a radius of the upwind sector that is 20 times greater than the average height of the roughness elements [13].
In this paper, upwind sectors with various dimensions and orientations were defined to perform a backscattering analysis and a correlation analysis with ground truth aerodynamic roughness lengths. The north direction was defined as the 0° orientation and the clockwise direction as the positive direction (see Figure 1(b)). Thirteen upwind sector radii were defined, in the range 1000 m to 7000 m, with steps of 500 m, and three upwind opening angles were defined: 30°, 45°, and 60°. Upwind sectors of the Beijing meteorological tower in the 25° orientation with an opening angle of 30° and radii of 1000 m, 3000 m, 5000 m, and 7000 m are illustrated in Figure 1(b). In our methodology, an upwind sector with specific dimensions was analyzed within the whole test image in 5° intervals, where the mean value of backscattering coefficient within the predefined upwind sector (illustrated using the green lines intersecting with the red circles in Figure 1(b)) was extracted and its correlation with the ground truth roughness length in the same direction was analyzed.
(4) The optimal dimension of the upwind sector was determined when the correlation between the backscattering coefficient within the sector and the ground truth roughness length in the same direction was at a maximum. Then the quantitative relationship model was set up between them, and the model was then used to map the aerodynamic roughness length in Beijing city using SAR data.
4. Results and Discussions
4.1. Backscattering Coefficient Analysis of the SAR Data
The backscattering coefficient within each upwind sector of varying dimensions and orientation, where the Beijing meteorological tower was taken as the central point, was analyzed. Figure 5 gives the results extracted from the ALOS PALSAR data on 21 June 2007, where we show sectors with a radius of 4500 m, opening angles of 30°, 45°, and 60°, respectively, and the orientation interval of the central angle of the upwind sector is 5°. Figure 5(a) is the ALOS PALSAR image on 21 June 2007, while Figures 5(b)–5(d) show the backscattering coefficient within the predefined upwind sector with opening angles of 30°, 45°, and 60°, respectively.
(a)
(b)
(c)
(d)
From Figure 5, it can be seen that there are great variations in the backscattering coefficient for various orientations of the Beijing meteorological tower. According to Figures 5(b), 5(c), and 5(d), the backscattering coefficient in the southwest (with an orientation of 180–270°) is higher because there are more tall and densely distributed buildings in this direction, whereas the backscattering coefficient in the northeast (0–90°) is lower because there are many green vegetation areas, and the buildings in this direction are sparsely distributed, as seen in Figure 5(a). Moreover, the backscattering coefficient extracted from the upwind sectors with an opening angle of 30° reveals more detailed information of the backscattering characteristics of the underlying surface at various orientations through comparing Figures 5(b), 5(c), and 5(d).
As mentioned above, the radius of the upwind sector also affects the aerodynamic roughness length calculation because there are different roughness elements contained in the sector for different radii. Thus the backscattering coefficients within the upwind sectors with radii changing from 1000 m to 7000 m, in steps of 500 m, were extracted and analyzed, where we took the Beijing meteorological tower as the center. Our results show that the backscattering coefficients extracted from various upwind sectors randomly change with radius because the urban areas are highly heterogeneous and the underlying roughness elements differ greatly when the radius is different. Thus in order to determine the optimal dimension of the upwind sector for the aerodynamic roughness calculation, the wind gradient data measured by the meteorological tower is needed for further analysis.
4.2. The Characteristics of the Aerodynamic Roughness Length Calculated Using the Wind Gradient Data
For the 106 day wind gradient data, experiments were firstly performed to evaluate whether they satisfied the neutral conditions or not according to the Richardson Number. And 86 groups of wind data that satisfied the neutral conditions requirement were used in the aerodynamic roughness length calculation via the fitting iteration method. Figure 6 illustrates the aerodynamic roughness length in each orientation, where the Beijing meteorological tower was taken as the center. According to Figure 6, the aerodynamic roughness length calculated using the wind gradient data ranged from 0 m to 6 m, where most were between 1 m and 3 m. It can be seen in Figure 6 that the value of to the southwest (180–240°) of the Beijing meteorological tower is larger because there are many buildings higher than 50 m in this direction, which is similar to the conclusions of Li et al. [17]. It can be seen in Figure 6 that the prevailing wind direction in Beijing city is from the northwest, because there are more data arising from 330–360° directions.
After calculation of from wind gradient data, we compared it with previous studies in the similar upwind direction. Then all the values of in the similar upwind direction calculated by us as well as those from literatures were plotted along with time. In Figure 7, our calculated value of was compared with the results of previous studies by Li et al. [17], Hu [20], Zhang and Chen [21], Gao et al. [22], and Ju [23]. Hu, Gao, and Zhang and Chen divided the range 0–360° into eight subranges, and Figure 7 shows within the 315–360° directions. Li divided 0–360° into four sections, and Figure 7 shows in the northwest direction (270–360°). For comparison in Figure 7 we show our calculated within the 330–360° direction. From Figure 7, it can be seen that the value of calculated by us is in good accordance with other studies, except that of Li, because the direction range considered by Li is much wider than ours, and the subsequent roughness elements in the various orientations are different. According to Figure 7, the aerodynamic roughness length to the northwest direction of the Beijing meteorological tower increased from 1985 to 2010, in accordance with the actual situations of Beijing city. And if a linear fitting is performed for all the values of except that of Li et al., the coefficient of determination is 0.76, and the root mean squared error is 0.35 m.
4.3. Correlations Analysis of SAR Backscattering and the Ground Truth Aerodynamic Roughness Length
In order to analyze the effect of the upwind sector dimension on the roughness length calculation, the following correlation analysis between the backscattering coefficient extracted from SAR data and the ground truth roughness length calculated using wind gradient data was carried out. The ground truth roughness length calculated using meteorological data was defined as , where represents the wind direction of the valid aerodynamic roughness length value numbered as ranging from 1 to 86. Then for the 22 scenes of ALOS PALSAR images, a series of upwind sectors with different dimensions were defined according to the acquired date and wind direction of the ground truth aerodynamic roughness length value , and then the backscattering coefficients within the upwind sector were extracted. Here represents the central angle of the sector, which corresponds to the wind direction of the th ground truth aerodynamic roughness length, represents the radius of the sector, which ranges from 1000 m to 7000 m with steps of 500 m, and is the opening angle of the sector, with values of 30°, 45°, and 60°. The coefficient of determination is () between the ground truth aerodynamic roughness length and the backscattering coefficient extracted from upwind sectors with certain radii and opening angles were calculated, respectively. Figure 8 illustrates the coefficient of determination between the backscattering coefficient and the ground truth for upwind sector radii ranging from 1000 m to 7000 m and opening angles of 30°, 45°, and 60°.
According to Figure 8, the coefficient of determination increases with an increase in sector radius when the radius is less than 2500 m. When the sector radius is larger than 2500 m, the coefficient of determination decreases with an increase in sector radius. When the radius is less than 3500 m, the coefficient of determination for the 30° opening angle is generally larger than the other two angles, and when the radius is larger than 3500 m, the determination coefficient of the 60° opening angle is generally larger than the other two angles. Among all the data, the coefficient of determination has a maximum of 0.802 for the 2500 m radius and 30° opening angle, which means that the roughness elements within the upwind sector of this dimension affect the aerodynamic roughness most.
The dimension of the fetch area which can best represent the effect of the roughness elements in the upwind direction has been kept as a controversial question. In this paper we determined the optimal dimension of the upwind sector through statistical correlation analysis between the backscattering coefficient extracted from SAR data and the ground truth roughness length calculated using wind gradient data. According to our knowledge, there are still lots of works which need to be carried out before we can find the physical rule behind the statistical result, and any comparison and benchmark works are highly welcomed.
4.4. Quantitative Relationship and Aerodynamic Roughness Length Mapping
Then the total of 86 pairs of ground truth aerodynamic roughness length calculated using the wind gradient data and the corresponding backscattering coefficient extracted from the optimal upwind sectors of 2500 m radius and 30° opening angle were randomly divided into two groups. The first group contained 58 pairs of data and was used for quantitative relationship investigation and modeling. Figure 9 shows the scatter plot between the backscattering coefficient and .
According to Figure 9, the higher the backscattering coefficient is, the higher the aerodynamic roughness is. As such, an exponential formula can be used to describe the relationship between them; that is,For (3), the coefficient of determination , the estimate error is 0.53 m, and the value for no correlation hypothesis testing is 0, showing the correlation is significant.
And the other group contained 28 pairs of data and was used for validation of (3) and the retrieved aerodynamic roughness results. Figure 10 shows the ground truth and the retrieved using (3) from ALOS PALSAR data, and the root mean square error is 0.37 m.
For the total of 86 pairs of ground truth and the corresponding backscattering coefficient extracted from SAR data, we changed the selected 58 pairs of training data and the form of the scatter plot between and the backscattering coefficient kept unchanged, the coefficient of (3) changed, yet the exponential formula was always the best form to quantitatively describe the relationship between and the backscattering coefficient. And the exponential formula always passed the validation of the remaining data pairs. That is to say, the exponential relationship between and backscattering coefficient in the upwind sector always exists and the roughness length of other areas can be mapped using exponential formula similar to (3).
Figure 11 shows the retrieved aerodynamic roughness length in 30° and 330° direction using the ALOS PALSAR image on 27 April 2010. Figure 11(a) shows the ALOS PALSAR image on 27 April 2010, and Figures 11(b) and 11(c) show the corresponding retrieved aerodynamic roughness length in 30° and 330° direction, where the arrow represents the wind direction, that is, the central direction of the upwind sector. For Figure 11(b), an upwind sector with a central angle of 30°, a radius of 2500 m, and an opening angle of 30° (i.e., from the 15° to 45° direction) was firstly defined, the backscattering coefficient within these sectors was extracted, and the aerodynamic roughness in 30° direction was retrieved using (3). The aerodynamic roughness in 330° direction shown in Figure 11(c) can be mapped in the similar way.
(a) ALOS PALSAR image on April 27, 2010 (the same extent with Figure 1(a))
(b) The map in the 30° direction (with the opening angle 15–45°)
(c) The map in the of 330° direction (with the opening angle 315–345°)
In Figures 11(b) and 11(c), the retrieved aerodynamic roughness length falls within the range 0.7–3.4 m, which is consistent with the ground truth data. For regions with a red color in the map, the roughness length is larger because the upwind sector contains more buildings, while for blue pixels in the map the roughness length is smaller because the upwind sector contains more water or green lands. Here take positions A and B as two examples. Position A lies near the Olympic forest park. The aerodynamic roughness length at position A in the 30° direction is 0.86 m because the roughness elements contained in the upwind sector in this direction mainly include trees, grass, and water. And the aerodynamic roughness length at position A in the 330° direction is 1.42 m because there are some residential buildings in the northwest direction. Position B lies near the commercial district. Its northeast is the Dazhongsi international square, and its northwest is Chinese Academy of Agricultural Sciences. Thus the calculated roughness length of position B in the 30° direction is 2.77 m and 1.89 m in the 330° direction, where the difference of the underlying roughness elements induces different aerodynamic roughness values. Thus the mapped roughness length from SAR data can effectively reveal detailed information of urban areas.
4.5. Temporal Variation Analysis of Aerodynamic Roughness Length Using Long Term SAR Data
According to analysis based on ALOS PALSAR data, the optimal dimension of upwind sector and the form of quantitative model can be referred and the aerodynamic roughness estimation models can be established for other SAR sensors, which is very important for long term study using multisource SAR data. Total of 22 scenes of COSMO-SkyMed data and 220 days of wind gradient data acquired from the Beijing meteorological tower were used to establish the quantitative relation between the COSMO-SkyMed and . The relation model is , the determination coefficient , the root mean square error is 0.18 m, and the value for no correlation hypothesis testing is 0, showing the correlation is significant. The correlation between the COSMO-SkyMed data and is even higher than ALOS PALSAR, which also is a good proof for the proposed methodology based on ALOS PALSAR data. Then multitemporal ALOS PALSAR and COSMO-SkyMed from 2006 to 2011 were, respectively, used to retrieve using corresponding relation model, and the temporal variation of aerodynamic roughness length of Beijing city was mapped and analyzed.
The temporal variation of the retrieved aerodynamic roughness length in the northeast of Beijing city in the 60° direction (opening angle from the 45° to 75° direction, the radius of 2500 m) is shown in Figure 12. It shows that the retrieved slowly increased with time because the buildings in the upwind sector were continuously increasing. Besides, the retrieved takes on seasonal patterns, with higher values in summer and lower values in winter, because the leaves of plants in summer increased the drag to air flow. Yet the decrease of in the winter of 2009 is not obvious, and conversely the retrieved between December of 2009 and March of 2011 is continually slowly climbing, not considering the seasonal effect of plants. It can be explained combining the land use change evolvement of this region (in Figure 13). Before December of 2009 this region was made up of about 30% buildings and more than 60% plants (see Figure 13(a)). Since December of 2009 a lot of buildings within this region had been constructed massively and till December of 2010 buildings within this region had reached 80% (see Figures 13(b)–13(d)). Thus it can be seen that the time series of retrieved from long term SAR data can well reveal the evolvement of the land use types in urban areas.
(a) June 13, 2009, COSMO-SkyMed
(b) December 6, 2009, COSMO-SkyMed
(c) December 9, 2010, COSMO-SkyMed
(d) March 15, 2011, ALOS PALSAR
5. Conclusions
In this paper the aerodynamic roughness length estimation methodology from SAR data was introduced by using Beijing city as an example. Wind gradient data acquired by the Beijing meteorological tower was used to calculate the ground truth aerodynamic roughness length. Then, the optimal dimension of the upwind sector was determined via a correlation analysis between the backscattering coefficient extracted from the SAR data and the corresponding ground truth aerodynamic roughness length determined in the same time and the same direction. And then a quantitative relationship, used to retrieve aerodynamic roughness length from SAR data, was established and the spatial and temporal aerodynamic roughness length in Beijing city was mapped using long term ALOS PALSAR and COSMO-SkyMed data. Our results show the retrieved from SAR data could accurately reveal the evolvement of the land use types and the surface roughness elements. Thus the methodology proposed in this paper can provide more accurate input for precise meteorological and climate applications in urban areas. However the limitations of this methodology lie in two aspects. One is that wind gradient data will be needed for quantitative modeling for applications in other cities, and the other is that the fluctuation of the backscattering coefficient of SAR data acquired at different dates will affect the accuracy of the roughness estimation. Aimed at this problem, polarimetric SAR is being studied to seek more stable algorithms for aerodynamic roughness estimation by us because some polarimetric parameters such as double bounce components derived from Freeman and Yamaguchi polarimetric decomposition are more stable indicators to buildings, which are the dominant urban roughness elements.
Competing Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.
Acknowledgments
This research was financially supported by National Natural Science Foundation of China (Grant nos. 41671359, 61471358), the Knowledge Innovation Program of the Chinese Academy of Sciences (Grant no. KZCX2-EW-320), National High Technology Research and Development Program of China (863 program) (Grant no. 2011AA120403), the ALOS research program (PI1404), the TanDEM research program (OTHER6984), and the China Scholarship Council (201604910064). Thanks are due to Dr. Junxia Dou for providing information about the meteorological tower.