Estimating global land surface broadband thermal-infrared emissivity using advanced very high resolution radiometer optical data

Abstract

An algorithm for retrieving global eight-day 5 km broadband emissivity (BBE) from advanced very high resolution radiometer (AVHRR) visible and near-infrared data from 1981 through 1999 was presented. Land surface was divided into three types according to its normalized difference vegetation index (NDVI) values: bare soil, vegetated area, and transition zone. For each type, BBE at 8–13.5 µm was formulated as a nonlinear function of AVHRR reflectance for Channels 1 and 2. Given difficulties in validating coarse emissivity products with ground measurements, the algorithm was cross-validated by comparing retrieved BBE with BBE derived through different methods. Retrieved BBE was initially compared with BBE derived from moderate-resolution imaging spectroradiometer (MODIS) albedos. Respective absolute bias and root-mean-square error were less than 0.003 and 0.014 for bare soil, less than 0.002 and 0.011 for transition zones, and −0.002 and 0.005 for vegetated areas. Retrieved BBE was also compared with BBE obtained through the NDVI threshold method. The proposed algorithm was better than the NDVI threshold method, particularly for bare soil. Finally, retrieved BBE and BBE derived from MODIS data were consistent, as were the two BBE values.

Keywords:

• earth observation
• broadband emissivity
• global land surface satellite
• soil taxonomy
• advanced very high resolution radiometer
• moderate-resolution imaging spectroradiometer

1. Introduction

Land surface thermal-infrared broadband emissivity (BBE) is an essential parameter for simulating surface energy budgets (Cheng et al. 2013; Liang et al. 2010; Pequignot, Chedin, and Scott 2008; Sellers, Dickinson, and Randall 1997). Given the lack of effective BBE data-sets, constant emissivity assumption or simple parameterization schemes are used in current land surface and general circulation models (Bonan 2002; Jin and Liang 2006; Zhou et al. 2003). A number of researchers have generated several BBE data-sets on a global or regional scale using different methods. Wilber et al. established a global BBE (5–100 µm) with a 10'×10' spatial resolution for satellite retrieval of long-wave radiation by assigning constant emissivity values to International Geosphere–Biosphere Program (IGBP) surface types (Wilber, Kratz, and Gupta 1999). Ogawa, Schmugge, and Rokugawa (2008) mapped the monthly global BBE (8.0–13.5 µm) using the moderate-resolution imaging spectroradiometer (MODIS) narrowband emissivity product (approximately 5 km) and a North African BBE (8.0–13.5 µm) using the advanced spaceborne thermal emission and reflectance radiometer (ASTER) narrowband emissivity product (90 m) (Ogawa and Schmugge 2004; Ogawa, Schmugge, and Rokugawa 2008). Two drawbacks were observed. (1) Either the spatial or temporal resolution of these products is limited. BBE with finer spatial and temporal resolutions is useful for local-scale studies of surface energy balance and for validating coarse resolution data, thus facilitating better understanding of land – atmosphere interactions (Liang et al. 2010; Ogawa and Schmugge 2004). (2) These products cannot capture variations in natural land surface BBE. In Wilber et al. (1999) seasonal variation was neglected, particularly for vegetated surfaces, because each IGBP surface type had only one emissivity over the entire growth stage. MODIS emissivity product cannot reflect vegetation cover effect on emissivity of vegetated surfaces (Wang and Liang 2009). A number of pioneer studies have demonstrated the importance of satellite-derived realistic BBE in improving results of global climate model simulation (Jin and Liang 2006; Zhou et al. 2003). Thus, BBE with high spatial and temporal resolutions for global land surface is urgently needed.

Converting satellite-derived narrowband emissivity products is a common method for deriving BBE (Ogawa and Schmugge 2004; Tang et al. 2011; Wang et al. 2005). However, retrieving land surface emissivity (LSE) from radiometric measurements of a thermal infrared sensor is inherently an ill-posed problem, to which solutions are likely to be unstable (Li, Strahler, and Friedl 1999; Liang 2001). Numerous temperature and emissivity inversion algorithms have been proposed to solve the aforementioned problem since the 1980s (Cheng et al. 2008, 2010a, 2011; Dash et al. 2002; Gillespie et al. 1998, 2011; Li and Becker 1993; Liang 2001, 2004; Liu and Xu 1998; Wan and Li 1997; Wang et al. 2008; Yu et al. 2009); however, only a few algorithms have been adapted to produce global LSE products, such as the temperature and emissivity separation (TES) algorithm for ASTER (Gillespie et al. 1998) and the day/night algorithm for MODIS (Wan and Li 1997). Cheng et al. examined potential methods for retrieving global land surface BBE with high spatial and temporal resolutions (Cheng and Liang 2012). ASTER narrowband emissivity product derived through the TES algorithm is inefficient because of its long revisit time of 16 days. MODIS narrowband emissivity product retrieved through the day/night algorithm lacks extensive validation and sufficient recognition from the remote sensing community, thus making it unsuitable for producing global land surface BBE with high spatial and temporal resolutions. Cheng et al. investigated the accuracy of remote sensing BBE at different spectral ranges in estimating net surface long-wave radiation and found that BBE can achieve optimal accuracy at a spectral range of 8.0–13.5 µm (Cheng et al. 2013). Based on soil spectra from the Johns Hopkins University Spectral Library, BBE at 8.0–13.5 µm and the seven corresponding MODIS narrowband black-sky albedos were calculated. A significant linear relationship between BBE and the seven narrowband black-sky albedos was explored. This relationship was initially verified for bare soil (Cheng and Liang 2012) using spatially and temporally matched ASTER emissivity product and MODIS narrowband albedo product, and then for homogeneous vegetated areas (Ren et al. 2012). This relationship was used to retrieve global eight-day 1 km land surface BBE using MODIS albedo products from 2000 through 2010. A comparison between derived BBE with that converted from the ASTER narrowband emissivity product indicated that these two BBE products agree with each other. Furthermore, BBE was validated through ground measurements obtained from several field campaigns in the USA and China. Absolute difference was 0.02.

This paper extends and adopts an algorithm for BBE retrieval from MODIS data to derive BBE from advanced very high resolution radiometer (AVHRR) visible and near-infrared (VNIR) data to produce global eight-day 5 km land surface BBE from 1981 through 1999. The remainder of this paper is structured as follows. Section 2 introduces data used. Section 3 describes the method of BBE retrieval from AVHRR VNIR data. Section 4 presents analysis of results. Section 5 discusses cross-validation results. Section 6 presents a brief conclusion.

2. Data

2.1. Global land surface satellite emissivity

Global land surface satellite (GLASS) emissivity is a BBE (8–13.5 µm) product derived from AVHRR VNIR data and MODIS albedos using our newly developed algorithms (Cheng and Liang 2012; Ren et al. 2012). GLASS emissivity has two parts: (1) the global eight-day 1 km land surface BBE retrieved from MODIS albedos from 2000 through 2010 and (2) the global eight-day 5 km land surface BBE retrieved from AVHRR VNIR reflectance from 1981 through 1999. In the algorithm used to generate GLASS BBE from MODIS albedos, land surface was classified into five types according to normalized difference vegetation index (NDVI) threshold values: water, snow or ice, bare soil (0< NDVI ≤0.156), vegetated area (NDVI >0.156), and transition zone (0.1< NDVI <0.2). Areas of overlapping bare soil and transition zone as well as of transition zone and vegetated area were noted. BBEs of water and snow or ice were set as 0.985 based on a combination of BBE calculated from the emissivity spectrum in the ASTER spectral library (http://spclib.jpl.nasa.gov) and the MODIS University of California Santa Barbara Emissivity Library (http://www.icess.ucsb.edu/modis/EMIS/html/em.html), and BBE simulated through radiative transfer models (Cheng et al. 2010b). BBEs of bare soil, vegetated areas, and transition zones were formulated as the linear function of seven MODIS narrowband black-sky albedos. When NDVI was lower than 0.1 or higher than 0.2, the formula for bare soil or vegetated areas was used to calculate individual BBE values. In areas of overlapping bare soil and transition zone (0.1< NDVI ≤0.156), BBE was taken as the average of values calculated using formulas for bare soil and transition zones. By contrast, BBE for areas of overlapping transition zone and vegetated area (0.156< NDVI <0.2) was taken as the average of values calculated using formulas for transition zones and vegetated areas. BBE derived from MODIS albedos was validated by field measurements in the US and China, and the absolute difference was 0.02 (Cheng and Liang 2012).

2.2. AVHRR reflectance data

AVHRR data used in this study were obtained from the AVH09 Surface Reflectance Product (Version 3) generated by the land long-term data record project, the goal of which is to produce a consistent long-term data-set from AVHRR and MODIS for use in global change and climate studies (Pedelty et al. 2007). Temporal resolution of AVH09 is one day. Processed AVHRR (AVH09) observations were grouped into separate scientific data-sets (SDSs) in a single hierarchical data format file, covering a 0.05° spatial resolution in a latitude/longitude climate modeling grid (CMG). SDSs included surface reflectance for Channels 1 (0.5–0.7 µm), 2 (0.7–1.0 µm), and 3 (3.55–3.93 µm) of AVHRR and top-of-atmosphere brightness temperature for Channels 3 (3.55–3.93 µm), 4 (10.3–11.3 µm), and 5 (11.5–12.5 µm). Reflectance and brightness temperature were atmospherically corrected, along with Rayleigh scattering, ozone, water vapor, and aerosol.

2.3. Auxiliary data

Soil taxonomy and MODIS land cover product (MOD12C1) were the required auxiliary data. The former was used to identify different soil orders, whereas the latter was used to identify different vegetation types. A soil taxonomy map was downloaded from the US Department of Agriculture (USDA) Natural Resources Conservation Service (NRCS) (http://soils.usda.gov/use/worldsoils/mapindex/order.html), which was based on a reclassification of the 1994 Food and Agriculture Organization – United Nations Educational, Scientific and Cultural Organization Soil Map of the World combined with a soil climate map. Spatial resolution of soil taxonomy is approximately 0.0333° at 5400 pixels×10,800 pixels.

Soil taxonomy employed is shown in Figure 1. Twelve soil orders can be observed in the map excluding shifting sand, rocky land, and ice or glacier. MOD12C1 is a yearly 0.05° spatial resolution CMG global land cover product with three different land cover classification schemes. IGBP global vegetation classification scheme was used to discriminate between vegetated and non-vegetated land covers, and to identify different vegetation types.

Figure 1. USDA NRCS soil taxonomy and selected study areas. Other classes, such as rocky land, shifting sand, and ice or glacier, are not included in the map.

3. Methodology

A flowchart for BBE retrieval from AVHRR VNIR data is presented in Figure 2. Both AVHRR VNIR data and BBE derived from MODIS albedos from year 2000 were acquired. AVHRR VNIR data and BBE were spatially and temporally matched. The eight-day 1 km BBE with a sinusoidal grid projection was reprojected to 0.05° CMG, and then, daily reflectance of AVH09 Channels 1 and 2 were combined at an eight-day step. Soil taxonomy was also reprojected to 0.05° CMG. High-quality clear-sky reflectance-BBE pairs were extracted from spatially and temporally matched data for each land surface type, which was identified by NDVI calculation from the combined eight-day AVHRR VNIR reflectance.

Figure 2. Flowchart of relationship established between MODIS BBE and AVHRR reflectance.

Similar to the algorithm used to derive BBE from MODIS albedos, land surface was divided into three types according to NDVI. NDVI threshold values used were based on previous studies (Momeni and Saradjian 2007; Sobrino et al. 2008). A bare soil pixel was identified as 0 <NDVI ≤0.2, a transition zone pixel as 0.145≤ NDVI ≤0.243, and a vegetated area pixel as NDVI ≥0.2. Finally, nonlinear formulas for bare soil, vegetated areas, and transition zones were established using nonlinear regression with extracted reflectance-BBE pairs.

Selected study areas (black rectangles in Figure 1) for bare soil and transition zones were the same. Relatively large homogeneous areas were selected for each soil order to obtain more representative samples for bare soil. Geographical locations of data used are presented in Table 1. Reflectance-BBE pairs extracted from data of days 49, 65, 81, 97, 105, 145, 241, and 321 in year 2000 were used to establish nonlinear relationship between BBE and VNIR reflectance. Reflectance-BBE pairs extracted from data of day 209 in year 2000 were used to test the established nonlinear relationship.

Table 1. Geographic locations for data used to develop the BBE retrieval algorithm for both bare soils and transition zones.

Soil order	Alfisols	Andisols	Aridisols	Entisols	Gelisols	Inceptisols
Geographical location
Latitude	39–42°	18.5–20.5°	35–40°	16–23°	31–36°	33–37°
Longitude	−86 to −80°	−102 to −97°	−120 to −110°	5–20°	87–96°	114–117°
Soil order	Mollisols	Oxisols	Ultisols	Vertisols	Spodosols	Histisols
Geographical location
Latitude	43–53°	−27.5 to −22°	32–36°	14–26°	48–51°	50–53°
Longitude	−108 to −98°	−55 to −50°	−85 to −80°	74–80°	−92.5 to 89.5°	−86.6 to −82.8°

Based on MOD12C1 data, 49 sites were selected from MODIS land subsets. The distribution of selected sites is shown in Figure 3. For each site, all data in 2000 were extracted from the 1°×1° square region centered on this site. Among these sites, 40 were used to establish nonlinear relationship between BBE and VNIR reflectance, and 9 were used to test the established relationship. Land cover for both developing and testing sites included forest, cropland, grassland, savanna, and shrub.

Figure 3. Distribution of sites used for developing and testing the algorithm for vegetated areas.

4. Analysis of results

4.1. Bare soil

For bare soil, 501,815 samples were obtained from the study areas. One formula was derived for all samples. However, this formula could not fully characterize the relationship between BBE and reflectance. Thus, two formulas were derived: one for individual vertisols and another for the remaining 11 soil orders. Each formula and the coefficient for each variable were significant at the 0.05 confidence level. The formulas are expressed as follows:

(1)

(2)

where ε_{BB_s1} and ε_{BB_s2} are BBEs for vertisols and the remaining 11 soil orders, respectively; and R₁ and R₂ are reflectance of AVHRR Channels l and 2, respectively. The BBE scatter plot predicted through Equation (1) versus BBE derived from MODIS albedos for vertisols and the histogram of bias are shown in Figure 4. Predicted BBE was lower than 0.97. Most points were distributed around the 1:1 line. However, points were highly dispersive when BBE was lower than 0.96. Thus, overall fit was poor, that is, correlation was only 0.529. Bias (as small as 4e-5) was centered on zero and distributed in a narrow band. Root-mean-square error (RMSE) was 0.004.

Figure 4. Scatter plot and difference histogram of BBE derived from MODIS albedos and BBE calculated through Equation (1).

The BBE scatter plot predicted through Equation (2) versus BBE derived from MODIS albedos for 11 soil orders and the histogram of bias are shown in Figure 5. Most points were distributed around the 1:1 line. Correlation was 0.822. Bias (−0.001) was centered on zero and distributed in a narrow band. RMSE was 0.012.

Figure 5. Scatter plot and difference histogram of BBE derived from MODIS albedos and BBE calculated through Equation (2).

4.2. Vegetation areas

A total of 92,691 samples were obtained to establish nonlinear relationship between BBE and reflectance for vegetated areas. BBE of vegetated areas was formulated as follows:

(3)

where ε_{BB_v} is BBE for vegetated areas. The scatter plot of BBE predicted through Equation (3) versus BBE derived from MODIS albedos and the histogram of bias are shown in Figure 6. BBE was clearly underestimated when it was higher than 0.98. Correlation was 0.588. Bias (as small as 2e-5) was centered on zero and distributed in a narrow band. RMSE was 0.007.

Figure 6. Scatter plot and difference histogram of BBE derived from MODIS albedos and BBE calculated through Equation (3).

4.3. Transition zones

A formula was derived for all soil orders by applying the same method used for bare soil. However, this formula did not work well for all soil orders, and thus, two formulas were derived: one for vertisols and another for the remaining 11 soil orders. Each formula and variable coefficient was significant at the 0.05 confidence level. The formulas are expressed as follows:

(4)

(5)

where ε_{BB_t1} and ε_{BB_t2} are BBEs for vertisols and the remaining 11 soil orders, respectively. The BBE scatter plot predicted through Equation (4) versus BBE derived from MODIS albedos for vertisols and the histogram of bias are shown in Figure 7. Predicted BBE was lower than 0.97. Points were dispersive when BBE was lower than 0.96. Correlation was 0.51. Bias (as small as 2e-5) was centered on zero and distributed in a narrow band. RMSE was 0.004. The BBE scatter plot predicted through Equation (5) versus BBE derived from MODIS albedos for the 11 soil orders and the histogram of bias are shown in Figure 8. Most points were distributed around the 1:1 line. Correlation was 0.867; bias was −0.002; and RMSE was 0.012.

Figure 7. Scatter plot and difference histogram of BBE derived from MODIS albedos and BBE calculated through Equation (4).

Figure 8. Scatter plot and difference histogram of BBE derived from MODIS albedos and BBE calculated through Equation (5).

5. Cross-validation

It is highly important to validate or test the algorithms developed for the retrieval of a certain particular biogeophysical parameter from satellite or aircraft data by using field measurements (Sobrino et al. 2006; Yu et al. 2012). However, emissivity validation is a complex process because the pixel is mixed and remains non-isothermal under most natural conditions, especially at moderate or coarse spatial resolutions. To our knowledge, field-measured emissivity at 5 km is unavailable. Thus, validating a 5-km emissivity product through field measurements is impractical. Three methods were used to test the algorithm applied in this study. First, retrieved BBE was compared with BBE retrieved from MODIS data, and then, with BBE derived through the NDVI threshold method (Sobrino et al. 2008), which is widely used to determine narrowband emissivity from AVHRR VNIR data. Finally, retrieved global BBE from year 2000 was compared with BBE derived from MODIS data to verify its consistency.

5.1. Comparison with MODIS BBE

The established formulas for bare soil, transition zones, and vegetated areas were examined based on extracted test data in Section 3. Test results for bare soil and transition zones are presented in Figures 9 and 10. For bare soil, absolute bias and RMSE were less than 0.003 and 0.014, respectively. For transition zones, absolute bias and RMSE were less than 0.002 and 0.011, respectively. Test result for vegetated areas is shown in Figure 11, with a bias and RMSE of −0.002 and 0.005, respectively. Retrieved BBE agreed with BBE derived from MODIS data, thus indicating that the established formulas can reflect the relationship between BBE and VNIR reflectance.

Figure 9. Difference histograms of BBE derived from MODIS data and BBE calculated through Equations (1) (left) and (2) (right).

Figure 10. Difference histograms of BBE derived from MODIS data and BBE calculated through Equations (4) (left) and (5) (right).

Figure 11. Difference histograms of BBE derived from MODIS data and BBE calculated through Equation (3).

5.2. Comparison with the NDVI threshold method

Based on Botswana's radiometer-measured data, which correspond with AVHRR channel, Van de Griend and Owe found that LSE has a good logarithmic relationship with NDVI (Van de Griend and Owe 1993). Through further analysis of the relationship between LSE and NDVI, Olioso noted that such dependence is closely related to soil emissivity, component effective emissivity of leaves, canopy structure, optical features of leaves, solar position, and proportion of soil penetrated by sunlight, although it is insensitive to observational geometric conditions (Olioso 1995). Valor and Caselles used a vegetation cover method to calculate LSE and applied this method to more complex mixed pixels to obtain satisfactory results (Valor and Caselles 1996). Sobrino et al. proposed an NDVI threshold method for AVHRR data, the mathematical expressions of which are presented as follows (Sobrino et al. 2008; Sobrino, Raissouni, and Li 2001):

(6)

(7)

where ε₄ and ε₅ are emissivities for AVHRR Channels 4 and 5; is fractional vegetation cover; NDVI_s is NDVI of soil; and NDVI_v is NDVI of vegetation. The pixel represents soil when NDVI < NDVI_s. The pixel is a mixture of vegetation and soil (i.e. a partial vegetation cover) when NDVI_s ≤ NDVI ≤ NDVI_v. The pixel represents vegetation when NDVI > NDVI_v. A linear function was derived from the ASTER spectral library to convert AVHRR narrowband emissivity to BBE at 8–13.5 µm. This function is presented in Equation (8). R² was 0.60 and RMSE was 0.018.

(8)

We derived Channel 4 emissivity from the same data used in Section 5.1 and converted it into BBE through Equation (8). Derived BBE was compared with BBE retrieved from MODIS data. Comparison results are shown in Figures 11 and 12. For vertisols, BBE retrieved through the NDVI threshold method was closer to BBE derived from MODIS data, with a bias and RMSE of −0.002 and 0.007, respectively. However, a high bias and RMSE were observed for the remaining 11 soil orders, with values as high as 0.028 and 0.035, respectively. We also retrieved BBEs for partially and fully vegetated areas and compared values with BBE derived from MODIS data. Bias and RMSE were less than 0.004 and 0.009, respectively. A comparison of Figures 9, 11, 12, and 13 reveals that BBE retrieved through the newly developed nonlinear formulas was better than BBE derived through the NDVI threshold method, especially for bare soil.

Figure 12. Difference histograms of BBE derived from MODIS albedos and BBE derived through the NDVI threshold method for vertisols (left) and for the remaining 11 soil orders (right).

Figure 13. Difference histograms of BBE derived from MODIS albedos and BBE derived through the NDVI threshold method for partially (left) and fully vegetated areas (right).

5.3. Consistency between AVHRR BBE and MODIS BBE

Verifying the consistency of multi-source data products is a prerequisite in generating products based on a long-time series of data, to which BBE product is no exception. Global land surface BBE was retrieved through nonlinear functions developed from AVHRR VNIR data for days 57, 137, 233, and 333 of year 2000, the values of which were compared with those derived from MODIS albedos. Difference maps are presented in Figure 14. The difference in the Sahara desert and in desert and semi-arid areas in Northwest China was larger than in other areas. Two AVHRR VNIR channels were obtained to calculate BBE, whereas seven VNIR channels were obtained for MODIS. The information provided by the two sensors was asymmetrical; thus, the difference in retrieved BBEs was unavoidable, albeit very small. Mean values were −0.0008, −0.0006, −0.0005, and 0.0006. RMSEs were 0.008, 0.008, 0.007, and 0.009.

Figure 14. Difference map of global land surface BBE derived from AVHRR VNIR data and MODIS albedos for days 57, 137, 233, and 333 of year 2000.

6. Conclusions

Satellite-derived land surface BBE with high spatial and temporal resolutions is important for improving studies on land surface energy balance. In the present study, the algorithm for retrieving global land surface BBE from MODIS albedos was adapted to retrieve BBE from AVHRR VNIR data. Land surface was divided into three types according to NDVI: bare soil, vegetated area, and transition zone. For each type, BBE at 8–13.5 µm was formulated as a nonlinear function of AVHRR reflectance for Channels 1 and 2. For bare soil, two formulas represented BBE: one for vertisols and another for the remaining 11 soil orders. R ²s of the fitted formulas were 0.529 and 0.828; biases were 4e-5 and −0.001; and RMSEs were 0.004 and 0.012 for vertisols and the remaining 11 soil orders, respectively. Two formulas were derived for transition zones: one for vertisols and another for the remaining 11 soil orders. R ²s of the fitted formulas were 0.51 and 0.867; biases were 2e-4 and −0.002; and RMSEs were 0.004 and 0.012 for vertisols and the remaining 11 soil orders, respectively. One formula was derived for vegetated areas. Corresponding fitted formula values were 0.588, 2e-5, and 0.007.

Given the difficulties in validating coarse emissivity products, we chose to validate BBE retrieved through the proposed algorithm by comparing it with BBE derived through other methods. Retrieved BBE agreed with BBE derived from MODIS albedos. Absolute bias and RMSE were less than 0.003 and 0.014 for bare soil, less than 0.002 and 0.011 for transition zones, and −0.002 and 0.005 for vegetated areas. Comparison results of retrieved BBE versus BBE derived through the NDVI threshold method indicate that the proposed algorithm is better than the NDVI threshold method, especially for bare soil. Consistency between retrieved BBE and BBE derived from MODIS albedos was also verified. Therefore, the two BBEs were consistent. The difference was very small, with mean values of −0.0008, −0.0006, −0.0005, and 0.0006; and RMSEs of 0.008, 0.008, 0.007, and 0.009 for days 57, 137, 233, and 333 of year 2000, respectively. The proposed algorithm was used to produce a GLASS BBE product from 1981 through 1999, which was released to the public in November 2012.

Acknowledgements

The MODIS and AVHRR data are obtained from https://wist.echo.nasa.gov/api/. The soil taxonomy is obtained from http://soils.usda.gov/use/worldsoils/mapindex/order.html. This work was supported in part by the National High Technology Research and Development Program of China via Grant 2009AA122100, by the National Natural Science Foundation of China via Grant 40901167 and 41201331 and the Fundamental Research Funds for the Central Universities.