Estimation of High Spatial-Resolution Clear-Sky Land Surface-Upwelling Longwave Radiation from VIIRS / S-NPP Data

Surface-upwelling longwave radiation (LWUP) is an important component of the surface radiation budget. Under the general framework of the hybrid method, the linear models and the multivariate adaptive regression spline (MARS) models are developed to estimate the 750 m instantaneous clear-sky LWUP from the top-of-atmosphere (TOA) radiance of the Visible Infrared Imaging Radiometer Suite (VIIRS) channels M14, M15, and M16. Comprehensive radiative transfer simulations are conducted to generate a huge amount of representative samples, from which the linear model and the MARS model are derived. The two models developed are validated by the field measurements collected from seven sites in the Surface Radiation Budget Network (SURFRAD). The bias and root-mean-square error (RMSE) of the linear models are −4.59 W/m2 and 16.15 W/m2, whereas those of the MARS models are−5.23 W/m2 and 16.38 W/m2, respectively. The linear models are preferable for the production of the operational LWUP product due to its higher computational efficiency and acceptable accuracy. The LWUP estimated by the linear models developed from VIIRS is compared to that retrieved from the Moderate Resolution Imaging Spectroradiometer (MODIS). They agree well with each other with bias and RMSE of −0.15 W/m2 and 25.24 W/m2 respectively. This is the first time that the hybrid method has been applied to globally estimate clear-sky LWUP from VIIRS data. The good performance of the developed hybrid method and consistency between VIIRS LWUP and MODIS LWUP indicate that the hybrid method is promising for producing the long-term high spatial resolution environmental data record (EDR) of LWUP.


Introduction
The surface radiation budget (SRB) is an important indicator in the study of climate formation and change and environmental prediction, which plays a key role in the global matter, energy cycles and interactions between the surface and the atmosphere system [1][2][3].SRB is dominated by longwave radiation in the night and during most of the year in the polar regions [4,5].Surface-upwelling longwave radiation (LWUP, 4.0-100 µm), the sum of thermal radiation emitted by the surface and reflected atmospheric downward longwave radiation, is the main cause of surface cooling in the clear night sky and is also an indirect indicator of surface temperature [5].An accurate estimate of LWUP is one of the prerequisites for obtaining accurate weather forecasts, climate simulations, and land-surface process simulations.
Generally, we can obtain the LWUP using three approaches: field measurement, satellite remote sensing and model prediction.LWUP can be accurately measured with field instruments.However, field networks are sparsely distributed globally.Furthermore, field measurement can only represent a limited area.The spatial resolution of model prediction is relatively coarse [6].Remote sensing can provide various kinds of products with global coverage and horizontal spatial continuity, but the temporal resolution is always limited.High spatial-resolution LWUP is an important diagnostic parameter for mesoscale land surface and atmosphere models [7] and can also serve as a bridge for validating coarse resolution data [8].The Visible Infrared Imaging Radiometer Suite (VIIRS) is one of the key instruments onboard the Suomi National Polar-orbiting Partnership (S-NPP) satellite system and was successfully launched in 2011.The spatial resolution of VIIRS thermal-infrared channels is 750 m.Ignoring the relatively low temporal resolution of VIIRS, it is the preferred data source for estimating high spatial-resolution LWUP.
Over the past few decades, substantial efforts have been devoted to estimating LWUP at the surface using remotely sensed observations from space-borne platforms.These methods can be primarily divided into two categories: the temperature-emissivity method [7][8][9][10] and the hybrid method [11][12][13].Current operational satellite land-surface temperature and emissivity products facilitate the estimation of LWUP [14 -17], but the large uncertainties in the land-surface temperature and emissivity products limit its accuracy [18,19].For example, the study of Wang et al. indicated that the accuracy of the temperature-emissivity method is much lower than that of the hybrid method [7].As shown in Figure 1, the weighting function of thermal infrared channels located in narrow bands that are semi-transparent to atmospheric gases and thus sensitive primarily to emission from the surface ("atmospheric windows") peaks at the surface and contains the surface-emission information, so LWUP can be derived from top-of-atmosphere (TOA) radiance or brightness temperature directly [13,20].The hybrid method links the thermal-infrared TOA radiances or brightness temperatures with LWUP through comprehensive radiative transfer modeling and statistical regression.It is physically based and at the same time has a high computational efficiency.Furthermore, it can bypass the problem of temperature and emissivity separation and achieve an acceptable accuracy in practice.The hybrid method has been successfully used to produce high spatial-resolution regional [7,11,21] and global [12] LWUP recently.
Remote Sens. 2018, 10, x FOR PEER REVIEW 2 of 16 cooling in the clear night sky and is also an indirect indicator of surface temperature [5].An accurate estimate of LWUP is one of the prerequisites for obtaining accurate weather forecasts, climate simulations, and land-surface process simulations.Generally, we can obtain the LWUP using three approaches: field measurement, satellite remote sensing and model prediction.LWUP can be accurately measured with field instruments.However, field networks are sparsely distributed globally.Furthermore, field measurement can only represent a limited area.The spatial resolution of model prediction is relatively coarse [6].Remote sensing can provide various kinds of products with global coverage and horizontal spatial continuity, but the temporal resolution is always limited.High spatial-resolution LWUP is an important diagnostic parameter for mesoscale land surface and atmosphere models [7] and can also serve as a bridge for validating coarse resolution data [8].The Visible Infrared Imaging Radiometer Suite (VIIRS) is one of the key instruments onboard the Suomi National Polar-orbiting Partnership (S-NPP) satellite system and was successfully launched in 2011.The spatial resolution of VIIRS thermal-infrared channels is 750 m.Ignoring the relatively low temporal resolution of VIIRS, it is the preferred data source for estimating high spatial-resolution LWUP.
Over the past few decades, substantial efforts have been devoted to estimating LWUP at the surface using remotely sensed observations from space-borne platforms.These methods can be primarily divided into two categories: the temperature-emissivity method [7][8][9][10] and the hybrid method [11][12][13].Current operational satellite land-surface temperature and emissivity products facilitate the estimation of LWUP [14-17], but the large uncertainties in the land-surface temperature and emissivity products limit its accuracy [18,19].For example, the study of Wang et al. indicated that the accuracy of the temperature-emissivity method is much lower than that of the hybrid method [7].As shown in Figure 1, the weighting function of thermal infrared channels located in narrow bands that are semi-transparent to atmospheric gases and thus sensitive primarily to emission from the surface ("atmospheric windows") peaks at the surface and contains the surface-emission information, so LWUP can be derived from top-of-atmosphere (TOA) radiance or brightness temperature directly [13,20].The hybrid method links the thermal-infrared TOA radiances or brightness temperatures with LWUP through comprehensive radiative transfer modeling and statistical regression.It is physically based and at the same time has a high computational efficiency.Furthermore, it can bypass the problem of temperature and emissivity separation and achieve an acceptable accuracy in practice.The hybrid method has been successfully used to produce high spatial-resolution regional [7,11,21] and global [12] LWUP recently.VIIRS was developed based on the heritage of Moderate Resolution Imaging Spectroradiometer (MODIS) instruments and has become a key bridge to ensure long-term continuity of the environmental data records (EDRs).VIIRS provides a large number of EDRs, including aerosol optical thickness [22], vegetation index [23], land-surface albedo [24], land-surface temperature [25], sea-surface temperature [26], etc.To our knowledge, an operational LWUP product from VIIRS is not available at regional and global scales.VIIRS do not provide operational products of surface broadband emissivity (BBE) and surface downward longwave radiation.Thus, it is difficult to calculate the LWUP with the temperature-emissivity method.In addition, no operational algorithm that can be used to retrieve LWUP from VIIRS has been reported in the literature.
We have developed a hybrid method for retrieving clear-sky LWUP from MODIS data and produced two years' global LWUP product recently [12].It is possible to produce long-term high spatial-resolution EDR of LWUP by combining MODIS and VIIRS.The first step is to develop a hybrid method to estimate the global 750-m instantaneous clear-sky LWUP from VIIRS data.The article is arranged as follows.The data including satellite data, field measurements, and atmospheric profiles are described in Section 2. The method and validation results are provided in Sections 3 and 4. A brief discussion and conclusion are given in Sections 5 and 6.

Visible Infrared Imaging Radiometer Suite (VIIRS) Data
VIIRS has been developed based on the heritage of legacy instruments, including AVHRR and MODIS, and extends and improves on them [27].The VIIRS has 22 spectral channels with wavelengths ranging from 0.41 µm to 12.5 µm, which can be used for environmental monitoring and numerical weather forecasting.The on-orbit verification and intensive calibration and validation using a ground target show that VIIRS is working very well [28,29].More than 20 environmental data records have been produced operationally from VIIRS data, including clouds, land-surface temperature and sea-surface temperature, vegetation index, aerosol optical thickness, active fire, snow/ice, surface albedo, etc.
The VIIRS data utilized in this study include the VIIRS sensor data records (SDR) and VIIRS Cloud Mask intermediate product (VCM).The VIIRS SDR contains the day-night band, imagery band, moderate resolution band and geolocation data.Three thermal-infrared channels, M14, M15 and M16, which are located in the "atmospheric windows" and are sensitive to the LWUP, are finally selected.These three channels have similar channel characteristics to MODIS channels 29, 31 and 32.The relative spectral responses of VIIRS channels M14, M15 and M16 as well as MODIS channels 29, 31 and 32 are displayed in Figure 1.The VIIRS VCM product is used to identify whether a pixel is clear-sky or cloudy [30].The pixels with confident clear flag are identified as clear-sky pixels.In addition, the longitude, latitude and the sensor view zenith angle data are also provided in the SDR.

Atmospheric Profiles
In this study, two years of the Atmospheric Infrared Sounder (AIRS) level 2 standard atmospheric profiles are collected to construct the atmospheric profile database.Launched aboard the Earth Observing System (EOS) second satellite Aqua, the AIRS instrument has 2378 infrared spectral channels covering 3.74-15.39µm with a high spectral resolution (λ/∆λ) of 1200 [31].The AIRS infrared band is very stable, and the offset is less than 10 mK/yr [32].The accuracy of temperature measurement of AIRS is better than 250 mK, and the absolute calibration accuracy of most bands can reach 100 mK.It is the most accurate and stable hyperspectral infrared detector up until now (http://daac.gsfc.nasa.gov/AIRS/).AIRS can be used to obtain the global atmospheric three-dimensional physical state (atmospheric temperature, water vapor, clouds, etc.) and the distribution of trace gases (ozone, carbon dioxide and methane, etc.) daily [33].The temperature profile of AIRS has 28 layers, and the corresponding atmospheric pressure varies from 1100 hPa to 0.1 hPa.Meanwhile, the AIRS water-vapor profile has 14 layers, and the corresponding atmospheric pressure ranges from 1100 hPa to 50 hPa.

Surface Radiation Budget Network (SURFRAD)
The Surface Radiation Budget Network (SURFRAD) network was established in 1993 to supply accurate, longstanding and persistent measurements of the surface radiation budget for climate change studies [34].SURFRAD provides quality-controlled field measurements including downward and upwelling solar irradiance and longwave infrared irradiation, along with other meteorological observations, such as wind speed, atmospheric pressure and relative humidity, etc. SURFRAD measurements are widely used for the validation of satellite-derived products and the details about the SURFRAD site and related instruments can be found in [34].Daily one or three-minute SURFRAD data are organized into daily ASCII text and can be freely downloaded from ftp://aftp.cmdl.noaa.gov/data/radiation/surfrad.In addition, the basic data are routinely sent to several archives including the World Radiation Data Center (WRDC), National Climatic Data Center (NCDC), and Baseline Surface Radiation Network (BSRN).The location, elevation and land cover of the seven SURFRAD sites are listed in Table 1.
Site upwelling and downward thermal infrared irradiance are measured using two precision infrared radiometers (PIR).The PIRs are sensitive to the spectral range from 3 µm to 50 µm.Three standard PIRs, which are annually calibrated by world-reputable organizations, are used to calibrate these two PIRs adopted in SURFRAD network.The spectral range of the measurements can be extended to 4-100 µm by calibration [35].The overall accuracy of PIR ground measurement is approximately ±9 W/m 2 [34], and is reported to be about ±5 W/m 2 recently [36].The VIIRS moderate resolution band (M-band) has a resolution of 750 m at nadir.The spatial matching issue and scale effect need be considered when validating satellite-derived LWUP using SURFRAD ground measurements.The footprint of SURFRAD PIRs is much smaller compared to that of the VIIRS.Fortunately, SURFRAD sites were selected at the locations where the surrounding land cover of the site was homogeneous.For example, Wang et al. [7] compared the brightness temperature of the ASTER pixel that contains the SURFRAD site to other neighboring pixels within 1 km × 1 km and 2 km × 2 km windows.They found that the discrepancies between the central pixel and surrounding pixels were less than 1 K in general, and the standard deviations of center pixel and surrounding pixels were less than 2 K under most conditions.Therefore, the SURFRAD measurements were used to validate the derived VIIRS LWUP directly.

Method
As shown in Figure 1, the channel characteristics of VIIRS channels M14, M15 and M16 are similar to those of MODIS channels 29, 31, and 32, which have been successfully used to derive LWUP using hybrid method at regional and global scales [7,12].In this section, we developed the hybrid method for VIIRS using TOA radiance of channels M14, M15 and M16 under the general framework of the hybrid method.First, a huge amount of representative samples are generated by extensive radiative transfer modeling; then the linear model and MARS model are established to predict LWUP using TOA radiance of channels M14, M15 and M16.

Radiative Transfer Modeling
The surface upwelling longwave radiation consists of two components: longwave radiation emitted by the surface and surface-reflected atmospheric downward longwave radiation, which can be written as: where the ε is the surface broadband emissivity.T s is the surface temperature.F l down is the downward longwave radiation.λ 1 and λ 2 are the spectral integration range (4-100 µm).
A representative simulation database is required to train the hybrid method, such as the linear model and the MARS model.The moderate resolution atmospheric transmission (MODTRAN) software [37], which is widely used by researchers all over the world, can be used to simulate the radiation transmission and interaction between the atmosphere and the land surface under various atmospheric and surface conditions.Allowing for the difference between land-surface temperature and atmospheric temperature, the global land surface is divided into three regions: the low-latitude region (30 [12].At the same time, the atmospheric profiles are also divided into these three categories according to the latitude.For each sub-region, we extracted the atmospheric profiles from AIRS Level 2 standard atmosphere product and constructed the atmospheric profile database.To avoid excessive computation and to alleviate the similarity of the profiles, a screening process was applied, and the criteria proposed in [12] adopted to measure the similarity of atmospheric profiles.In total, 2842, 35,487 and 41,724 atmospheric profiles of high-latitude, mid-latitude and low-latitude region were obtained.MODTRAN 5.2 was used to simulate the spectral downward longwave radiance, thermal path radiance and spectral transmittance for each atmospheric profile and sensor view zenith angle.During the simulation, the sensor view zenith angles were set from 0 • to 60 • with an interval of 15 • .We calculated the difference between surface temperature and the bottom layer temperature of atmospheric profiles for each region using two-year AIRS standard L2 product.The surface temperature was determined based on this difference.For example, the surface temperature of mid-latitude is equal to the bottom layer temperature plus a range of [−10, 15] K with a step of 5 K. Eighty-four representative spectral emissivity spectra including vegetation, soil, snow/ice, water selected from ASTER spectral library [38] and MODIS UCSB spectral library [39] and their combinations were used to characterize the land-surface conditions.Since there was no spectral emissivity value for wavelengths larger than 14 µm, the emissivity value was supposed to be the same as that at wavelength 14 µm in the following simulations.When the surface and atmospheric parameters were determined, the VIIRS TOA channel radiances are calculated using the following simplified equation: where L i is the TOA radiance for channel M14, M15 and M16 of VIIRS; L ↑λ and L ↓λ are path thermal radiance and spectral downward longwave radiance, respectively; f i (λ) is the spectral response function for channels M14, M15 and M16; λ 1 and λ 2 are the spectral range of VIIRS TIR channel.The flowchart of the hybrid method is displayed in Figure 2.

Linear Model
A linear model with the following expression was developed using the generated samples in Section 3.1: where 0 a , 1 a , 2 a and 3 a are coefficients; and M14, M15 and M16 are TOA radiance for moderate resolution bands 14, 15 and 16 of the VIIRS.The linear model is fitted for each sub-region and view zenith angle.

The Multivariate Adaptive Regression Spline (MARS) Model
To probe the non-linear relationship between TOA spectral radiance and LWUP, we also use the multivariate adaptive regression spline (MARS) to model the non-linear relationship between TOA spectral radiance and LWUP with the same samples as those used by the linear model.MARS was proposed by Jerome H. Friedman in 1991 [40].MARS is a highly generalized and highly specialized regression method for high-dimensional data.The regression method takes the tensor product of spline functions as the basis functions, while the determination of the basis functions and their number are automatically completed by the data without manual selection.In the multi-dimensional case, how to divide the space has become a critical problem, but the MARS model can solve this problem well.The MARS model is defined as follows:

Linear Model
A linear model with the following expression was developed using the generated samples in Section 3.1: where a 0 , a 1 , a 2 and a 3 are coefficients; and M14, M15 and M16 are TOA radiance for moderate resolution bands 14, 15 and 16 of the VIIRS.The linear model is fitted for each sub-region and view zenith angle.

The Multivariate Adaptive Regression Spline (MARS) Model
To probe the non-linear relationship between TOA spectral radiance and LWUP, we also use the multivariate adaptive regression spline (MARS) to model the non-linear relationship between TOA spectral radiance and LWUP with the same samples as those used by the linear model.MARS was proposed by Jerome H. Friedman in 1991 [40].MARS is a highly generalized and highly specialized regression method for high-dimensional data.The regression method takes the tensor product of spline functions as the basis functions, while the determination of the basis functions and their number are automatically completed by the data without manual selection.In the multi-dimensional case, how to divide the space has become a critical problem, but the MARS model can solve this problem well.The MARS model is defined as follows: where B m (x) is the mth basis function; a 0 is coefficient; a m is the coefficient of the mth basis function; M is the number of the basis function; K m is the number of knots in the mth basis function; S km is 1 or −1, and indicates the spline function on the right or left side;v(k, m) labels the independent variables; and t km is a knot location.The basis function of MARS is a single spline function, or is the result of the interaction of multiple spline functions.The spline functions on the right and left sides are defined as follows: where t is the position of the node; x − t km and t km − x are used to describe spline functions on the right and left regions when t is given; q is the power (>0) to which the splines are raised; and "+" takes 0 for negative values.MARS is built in two phases: the forward-basis function selection and backward-pruning process.In the first stage, MARS begins with just the intercept term and iterations are conducted to add pairs of basis functions, which can minimize the training error to the utmost extent.In order to avoid overfitting, a backward-pruning process is needed.In the backward-pruning phase, the basis function, which contributes the least to the reduction of training error, is deleted at a time.The model with the lowest GCV (generalized cross-validation) is finally selected.The GCV, as an estimator for effective mean residual error, can achieve a balance between goodness of fit and model complexity.The GCV is calculated as follows: where N is the number of samples in the training data; f (x i ) is the estimation of y i ; and enp is the effective number of parameters: where c is the penalty parameter; and k is the number of non-constant basis functions.
The training of MARS is based on the ARESLab package of the MATLAB platform, in which the MARS is implemented according to Friedman's original papers [40] and all parameters in the ARESLab package are automatically determined.The source code of the ARESLab toolbox can be downloaded from the following website: www.cs.rtu.lv/jekabsons/regression.html.

Training Results of the Linear Model and MARS Model
The linear model can account for more than 97.7%, 98.5% and 99.1% of the variation of the LWUP in the simulation database for the low-latitude, middle-latitude and high-latitude regions, respectively.The bias is zero, and the RMSE ranges from 5.27 to 13.02 W/m 2 .The RMSE in the low-latitude region is larger than that in the middle-latitude and high-latitude region.In addition, the RMSE for the larger view zenith angle is greater than that with the smaller view zenith angle.Details about fitting the results of the linear models can be found in Table 2.
The training results of the MARS model are also displayed in Table 2.The MARS model can account for more than 98.2%, 98.7% and 99.2% of the variation in the simulation database for the low-latitude, middle-latitude and high-latitude regions, respectively.The biases of all MARS models are zero, and the RMSEs range from 8.53 to 10.13 W/m 2 , 6.64 to 8.24 W/m 2 , and 4.8 to 6.47 W/m 2 for the low-latitude, middle-latitude and high-latitude regions, respectively.The RMSEs of the MARS models are slightly less than those of the linear models.The field measurements collected from the SURFRAD network were used to validate the linear models developed.When the view zenith angle was not equal to 0 • , 15 • , 30 • , 45 • , 60 • , LWUP was linearly interpolated from LWUPs predicted by the linear model with an adjacent view zenith angle.View zenith angle exceeding 60 • was not considered.The cloud-mask information was extracted from the VIIRS Cloud Mask intermediate product.Clear sky was identified when 3 × 3 neighboring pixels of the site are clear to ensure that the field measurements at the site were not affected by clouds and cloud shadows.In total, 2901 validation samples were finally obtained.
Figure 3 shows the comparison between the estimated LWUP and the ground measurements.We can find that the average bias and RMSE are −4.59 and 16.15 W/m 2 , respectively, at seven SURFRAD sites.For further analysis, the validation samples were divided into two groups: daytime and night time according to the local solar time.The bias ranged from −16.94 to 7.84 W/m 2 , and RMSE ranges from 13.84 to 30.43 W/m 2 in the daytime.The average bias and RMSE were −1.40 W/m 2 and 21.57W/m 2 .If Desert_Rock_NVt is not considered, the average bias and RMSE were −1.12 W/m 2 and 21.71 W/m 2 .During the night time, the bias ranged from −20.48 W/m 2 to 2.97 W/m 2 , and RMSE ranged from 8.76 W/m 2 to 21.22 W/m 2 .The average bias and RMSE were −8.3 W/m 2 and 13.46 W/m 2 .If Desert_Rock_NVt is not considered, the average bias and RMSE were −5.33 W/m 2 and 10.73 W/m 2 .The LWUP estimated at night was less divergent than that in the daytime.

MARS Models
The LWUP estimated by MARS models was validated with the same ground measurements.As shown in Figure 4, we can find that the average bias and RMSE were −5.23 and 16.38 W/m 2 .The bias ranged from −16.15 to 8.15 W/m 2 , and RMSE ranged from 12.12 to 29.96 W/m 2 in the daytime.The

MARS Models
The LWUP estimated by MARS models was validated with the same ground measurements.As shown in Figure 4, we can find that the average bias and RMSE were −5.The average bias and RMSE were −8.93 W/m 2 and 13.8 W/m 2 .If Desert_Rock_NVt is not considered, the average bias and RMSE were −5.93 W/m 2 and 11.03 W/m 2 .The distribution of the points in Figure 4 is similar to that in Figure 3.
Regarding the satellite-derived surface radiative fluxes, an accuracy of approximately ±20 W/m 2 for instantaneous footprint values is required by the hydrological, meteorological, and agricultural research communities [41].According to the validation results in this study, both the linear and MARS models can meet this requirement.Overall, the linear models are slightly better than the MARS models, although MARS can model the non-linearity between LWUP and TOA spectral radiance well during the training stage.Two reasons may account for this phenomenon.First, the zenith angles of the validation samples are all less than 40 • , and the differences between the linear models and the MARS models are slight when the zenith angle is less than 45 • , as shown in Table 2. Second, various kinds of satellite measurement errors are likely to be magnified by the non-linear model such as MARS.Due to the higher computational efficiency of the linear models, it is more adaptable to produce operational LWUP product.
A few studies have been devoted to estimating high spatial-resolution LWUP from MODIS and VIIRS.For example, Wang et al. [7] developed the linear models for estimating North American LWUP using MODIS data.The average bias and RMSE over SURFRAD site were −10.97 W/m 2 and 18.35 W/m 2 .Jiao et al. [11] developed the neural network models for estimating LWUP from MODIS and VIIRS data over the Tibet Plateau.The average bias and RMSE of the validation results for MODIS were 11.24 W/m 2 and 26.78 W/m 2 .The accuracy of LWUP from VIIRS was not validated.Cheng and Liang developed the linear models for estimating LWUP from MODIS data at global scale.The bias and RMSE of the validation over SURFRAD site were −4.49W/m 2 and 13.47 W/m 2 .Compared to these studies, the accuracies of the newly developed hybrid method are comparable or superior to those of the related published references.
average bias and RMSE were −1.27 W/m 2 and 19.79 W/m 2 .If Desert_Rock_NVt is not considered, the average bias and RMSE were −1.13 W/m 2 and 21.64 W/m 2 .During the nighttime, the bias ranged from −20.46 W/m 2 to 2.25 W/m 2 , and RMSE ranged from 8.66 W/m 2 to 15.26 W/m 2 .The average bias and RMSE were −8.93 W/m 2 and 13.8 W/m 2 .If Desert_Rock_NVt is not considered, the average bias and RMSE were −5.93 W/m 2 and 11.03 W/m 2 .The distribution of the points in Figure 4 is similar to that in Figure 3.
Regarding the satellite-derived surface radiative fluxes, an accuracy of approximately ±20 W/m 2 for instantaneous footprint values is required by the hydrological, meteorological, and agricultural research communities [41].According to the validation results in this study, both the linear and MARS models can meet this requirement.Overall, the linear models are slightly better than the MARS models, although MARS can model the non-linearity between LWUP and TOA spectral radiance well during the training stage.Two reasons may account for this phenomenon.First, the zenith angles of the validation samples are all less than 40°, and the differences between the linear models and the MARS models are slight when the zenith angle is less than 45°, as shown in Table 2. Second, various kinds of satellite measurement errors are likely to be magnified by the non-linear model such as MARS.Due to the higher computational efficiency of the linear models, it is more adaptable to produce operational LWUP product.
A few studies have been devoted to estimating high spatial-resolution LWUP from MODIS and VIIRS.For example, Wang et al. [7] developed the linear models for estimating North American LWUP using MODIS data.The average bias and RMSE over SURFRAD site were −10.97 W/m 2 and 18.35 W/m 2 .Jiao et al. [11] developed the neural network models for estimating LWUP from MODIS and VIIRS data over the Tibet Plateau.The average bias and RMSE of the validation results for MODIS were 11.24 W/m 2 and 26.78 W/m 2 .The accuracy of LWUP from VIIRS was not validated.Cheng and Liang developed the linear models for estimating LWUP from MODIS data at global scale.The bias and RMSE of the validation over SURFRAD site were −4.49W/m 2 and 13.47 W/m 2 .Compared to these studies, the accuracies of the newly developed hybrid method are comparable or superior to those of the related published references.

Comparison between Moderate Resolution Imaging Spectroradiometer (MODIS) and VIIRS Surface-Upwelling Longwave Radiation (LWUP)
Remote-sensing instruments such as MODIS and VIIRS are critical for providing continual and reliable measurements of the atmosphere, ocean and land-surface variables at the global scale.Furthermore, VIIRS is developed based on the heritage of MODIS instruments and has become a key bridge to ensure long-term continuity of the climate data records.LWUP estimated from VIIRS is compared with that retrieved from MODIS using the linear models of Cheng et al. [12] to check the consistency between the two instruments.The overpass time of the MODIS granule image is at 20:40 (UTC), 17 August 2014 and the VIIRS granule image is acquired from 20:46 to 20:51 (UTC) on the same day.The time difference between the two images is less than 10 min.Therefore, it is supposed that there are no significant differences in atmospheric conditions and land-surface properties within the 10-min period.In addition, the corresponding MOD35 and VIIRS Cloud Mask intermediate products are used to eliminate the influence of clouds.The linear hybrid method proposed by Cheng et al. [12] is used to calculate the LWUP from MODIS data, which has been validated by three measurement networks with a bias and RMSE of −0.31 W/m 2 and 19.92 W/m 2 in total.
The VIIRS image was aggregated to 1 km with bi-linear resampling method to match the spatial resolution of MODIS.The LWUP derived from both MODIS and VIIRS data were compared without considering the pixels that covered by invalid data or were contaminated by cloud.The spatial distributions of LWUP for MODIS and VIIRS are displayed in Figure 5, from which we can find that the LWUP distribution pattern of VIIRS is similar to that of the MODIS.In addition, the density plot of the LWUP is shown in Figure 6.Most of the LWUP values are concentrated between 400 and 800 W/m 2 with R 2 of 0.88.The bias and RMSE are −0.15W/m 2 and 25.24 W/m 2 , respectively.Thus, the LWUP of VIIRS is consistent with that of MODIS.This result is consistent with the study of Jiao et al. [11].They compared the LWUP retrieved from MODIS and VIIRS at Tibet Plateau.The R 2 , bias and RMES were 0.52 W/m 2 , 2.87 W/m 2 and 26.02 W/m 2 , respectively.

Comparison between Moderate Resolution Imaging Spectroradiometer (MODIS) and VIIRS Surface-Upwelling Longwave Radiation (LWUP)
Remote-sensing instruments such as MODIS and VIIRS are critical for providing continual and reliable measurements of the atmosphere, ocean and land-surface variables at the global scale.Furthermore, VIIRS is developed based on the heritage of MODIS instruments and has become a key bridge to ensure long-term continuity of the climate data records.LWUP estimated from VIIRS is compared with that retrieved from MODIS using the linear models of Cheng et al. [12] to check the consistency between the two instruments.The overpass time of the MODIS granule image is at 20:40 (UTC), 17 August 2014 and the VIIRS granule image is acquired from 20:46 to 20:51 (UTC) on the same day.The time difference between the two images is less than 10 min.Therefore, it is supposed that there are no significant differences in atmospheric conditions and land-surface properties within the 10-min period.In addition, the corresponding MOD35 and VIIRS Cloud Mask intermediate products are used to eliminate the influence of clouds.The linear hybrid method proposed by Cheng et al. [12] is used to calculate the LWUP from MODIS data, which has been validated by three measurement networks with a bias and RMSE of −0.31 W/m 2 and 19.92 W/m 2 in total.
The VIIRS image was aggregated to 1 km with bi-linear resampling method to match the spatial resolution of MODIS.The LWUP derived from both MODIS and VIIRS data were compared without considering the pixels that covered by invalid data or were contaminated by cloud.The spatial distributions of LWUP for MODIS and VIIRS are displayed in Figure 5, from which we can find that the LWUP distribution pattern of VIIRS is similar to that of the MODIS.In addition, the density plot of the LWUP is shown in Figure 6.Most of the LWUP values are concentrated between 400 and 800 W/m 2 with R 2 of 0.88.The bias and RMSE are −0.15W/m 2 and 25.24 W/m 2 , respectively.Thus, the LWUP of VIIRS is consistent with that of MODIS.This result is consistent with the study of Jiao et al. [11].
They compared the LWUP retrieved from MODIS and VIIRS at Tibet Plateau.The R 2 , bias and RMES were 0.52 W/m 2 , 2.87 W/m 2 and 26.02 W/m 2 , respectively.

Cloud Effect
In Section 4, the established hybrid methods were validated by the field measurements.Generally, most of the estimated LWUPs were closer to the one-to-one line.The points were less divergent in the night time than in the daytime.This may be because that the land surface is more homogeneous in the night time than in the daytime [42].Compared to other sites, obvious underestimation is found in the desert site.Wang and Liang [7] thought air traffic out of Los Angeles

Cloud Effect
In Section 4, the established hybrid methods were validated by the field measurements.Generally, most of the estimated LWUPs were closer to the one-to-one line.The points were less divergent in the night time than in the daytime.This may be because that the land surface is more homogeneous in the night time than in the daytime [42].Compared to other sites, obvious

Cloud Effect
In Section 4, the established hybrid methods were validated by the field measurements.Generally, most of the estimated LWUPs were closer to the one-to-one line.The points were less divergent in the night time than in the daytime.This may be because that the land surface is more homogeneous in the night time than in the daytime [42].Compared to other sites, obvious underestimation is found in the desert site.Wang and Liang [7] thought air traffic out of Los Angeles produces many cirrus clouds over this site, and cloud contamination may be a significant source of error at this site.Thus, one of the reasons may be that the VIIRS Cloud Mask intermediate product cannot identify all kinds of cloud types, especially the cirrus cloud.
In order to test this assumption, we downloaded three years (2014-2017) field measurements from Atmospheric Radiation Measurement (ARM) program Southern Great Plains (SGP) site C1 (https://www.arm.gov).The downloaded data include the LWUP, cloud-base height and other auxiliary data.The cloud-base height information is derived from the Micropulse Lidar (MPL).MPL is a ground-based remote-sensing instrument that can be used to measure the altitude of clouds by transmitting pulses of light and using the receiver to detect the light scattered back by clouds and precipitation.From the time delay between each outgoing pulse and the backscattered signal, the distance to the scatterer is inferred [43].The VIIRS pixels are, first, identified by the VIIRS Cloud Mask intermediate product, and then the extracted clear-sky pixels are further screened by the field measurements.The pixels with more than one detected cloud base are considered as cloudy, while those with no significant backscatter detected are regarded as true clear-sky pixels.The developed liner models are used to retrieve LWUP.The validation results of LWUP at site C1 are displayed in

Broadband Emissivity (BBE) Effect
The obvious underestimation of LWUP at Desert_Rock_NV may also be related to BBE of the surface.Taking the developed linear model for the mid-latitude region at nadir view as an example, we investigated the effects of surface BBE on the accuracy of LWUP estimation.The bias and RMSE of the fitting liner model at nadir are zero and 6.94 W/m 2 .We calculated the average bias and RMSE for each emissivity spectrum as well as the corresponding BBE.The relationship between BBE versus bias is shown in Figure 8; 78% samples have a negative bias when their BBE is less than 0.966 and 88% BBE has a positive bias when BBE is equal to or larger than 0.966.It is very likely to produce negative bias at Desert_Rock_NV from the point of the developed linear models, because its BBE is less than 0.966.We have also investigated the relationship between the fitting residual and landsurface temperature (LST) with the same data, and no significant overestimated or underestimated trend is found.Thus, cloud and surface BBE are two primary factors that affect the accuracy of the

Broadband Emissivity (BBE) Effect
The obvious underestimation of LWUP at Desert_Rock_NV may also be related to BBE of the surface.Taking the developed linear model for the mid-latitude region at nadir view as an example, we investigated the effects of surface BBE on the accuracy of LWUP estimation.The bias and RMSE of the fitting liner model at nadir are zero and 6.94 W/m 2 .We calculated the average bias and RMSE for each emissivity spectrum as well as the corresponding BBE.The relationship between BBE versus bias is shown in Figure 8; 78% samples have a negative bias when their BBE is less than 0.966 and 88% BBE has a positive bias when BBE is equal to or larger than 0.966.It is very likely to produce negative bias at Desert_Rock_NV from the point of the developed linear models, because its BBE is less than 0.966.We have also investigated the relationship between the fitting residual and land-surface temperature (LST) with the same data, and no significant overestimated or underestimated trend is found.Thus, cloud and surface BBE are two primary factors that affect the accuracy of the LWUP estimate.

Conclusions
LWUP is an important component of the land-surface radiation budget.We developed two hybrid methods, namely, the linear model and MARS model, to estimate the 750-m instantaneous clear-sky LWUP from the VIIRS TOA channel radiances of M14, M15, and M16 at the global scale.
To consider the difference between the land-surface temperature and air temperature, the global land surface was divided into 3 regions based on latitude.Extensive radiative transfer modeling was conducted to produce a huge amount of representative samples for each sub region.The linear models and the MARS models were established at 0°, 15°, 30°, 45° and 60° viewing zenith angles for each sub region.According to the statistical results, the linear models can account for more than 97.7% of the variation in the simulation database, the bias is zero, and the RMSEs range from 5.27 to 13.02 W/m 2 ; the MARS models can account for more than 98.2% of the variation in the simulation database, the bias is zero, and the RMSEs range from 4.8 to 10.13 W/m 2 .Then, the two models were validated using three years (2014-2017) of ground measurements collected from seven SURFRAD sites.The average bias and RMSE of the linear models were −4.59 W/m 2 and 16.15 W/m 2 , whereas the average bias and RMSE of the MARS models were −5.23 W/m 2 and 16.38 W/m 2 .The difference between the linear models and the MARS models was not significant.The linear models have higher computational efficiency and are easy to implement, so it is a good choice for producing the operational LWUP product.The LWUP retrieved from VIIRS by the developed linear models was compared to that retrieved from MODIS by the previous linear models developed.Their spatial distribution pattern agreed well with each other; the bias and RMSE were −0.15 W/m 2 and 25.24 W/m 2 .This is the first time that the hybrid method has been applied to global estimates of clear-sky LWUP from VIIRS data.Limited validation results indicate that the accuracy of the hybrid method can meet the accuracy requirement of the hydrological, meteorological, and agricultural research communities.The good performance of the developed hybrid method and consistency between VIIRS LWUP and MODIS LWUP indicate that the hybrid method shows promise for producing a long-term high spatial-resolution environmental data record (EDR) of LWUP.The field measurements used for validation are collected from the SURFRAD network, which is located in the middle-latitude region.More validation data from low-latitude and high-latitude regions with other land-cover types need to be collected for further validation in the future.Also, more extensive comparison should be conducted before the generation of LWUP EDR using MODIS and VIIRS data.
be collected for further validation in the future.Also, more extensive comparison should be conducted before the generation of LWUP EDR using MODIS and VIIRS data.

Figure 1 .
Figure 1.Relative spectral responses of Visible Infrared Imaging Radiometer Suite (VIIRS) channels M14, M15 and M16, and Moderate Resolution Imaging Spectroradiometer (MODIS) channels 29, 31 and 32.The gray line represents the transmittance of the 1976 U.S. Standard Atmosphere.

Figure 1 .
Figure 1.Relative spectral responses of Visible Infrared Imaging Radiometer Suite (VIIRS) channels M14, M15 and M16, and Moderate Resolution Imaging Spectroradiometer (MODIS) channels 29, 31 and 32.The gray line represents the transmittance of the 1976 U.S. Standard Atmosphere.

Figure 2 .
Figure 2. Flowchart for developing the hybrid method.
the th m basis function; 0 a is coefficient; m a is the coefficient of the th m basis function; M is the number of the basis function; m K is the number of knots in the th m basis function; km S is 1 or −1, and indicates the spline function on the right or left side; v( , ) k m labels the independent variables; and km t is a knot location.The basis function of MARS is a single spline function, or is the result of the interaction of multiple spline functions.The spline functions on the right and left sides are defined as follows:

Figure 2 .
Figure 2. Flowchart for developing the hybrid method.
23 and 16.38 W/m 2 .The bias ranged from −16.15 to 8.15 W/m 2 , and RMSE ranged from 12.12 to 29.96 W/m 2 in the daytime.The average bias and RMSE were −1.27 W/m 2 and 19.79 W/m 2 .If Desert_Rock_NVt is not considered, the average bias and RMSE were −1.13 W/m 2 and 21.64 W/m 2 .During the nighttime, the bias ranged from −20.46 W/m 2 to 2.25 W/m 2 , and RMSE ranged from 8.66 W/m 2 to 15.26 W/m 2 .

Figure 6 .
Figure 6.The scatterplot of retrieved LWUP from VIIRS and that retrieved from MODIS.

Figure 5 .
Figure 5. Distribution of LWUPs derived from VIIRS and MODIS images.

Figure 6 .
Figure 6.The scatterplot of retrieved LWUP from VIIRS and that retrieved from MODIS.

Figure 6 .
Figure 6.The scatterplot of retrieved LWUP from VIIRS and that retrieved from MODIS.

Figure 7 .
The bias and RMSE of cloudy pixels are −6.65 W/m 2 and 16.04 W/m 2 , while bias and RMSE of clear-sky pixels are −2.94W/m 2 and 15.62 W/m 2 .The VIIRS Cloud Mask intermediate product actually cannot identify all kinds of cloud types, and undetected cloud will cause the underestimation of LWUP.Remote Sens. 2018, 10, x FOR PEER REVIEW 13 of 16produces many cirrus clouds over this site, and cloud contamination may be a significant source of error at this site.Thus, one of the reasons may be that the VIIRS Cloud Mask intermediate product cannot identify all kinds of cloud types, especially the cirrus cloud.In order to test this assumption, we downloaded three years (2014-2017) field measurements from Atmospheric Radiation Measurement (ARM) program Southern Great Plains (SGP) site C1 (https://www.arm.gov).The downloaded data include the LWUP, cloud-base height and other auxiliary data.The cloud-base height information is derived from the Micropulse Lidar (MPL).MPL is a ground-based remote-sensing instrument that can be used to measure the altitude of clouds by transmitting pulses of light and using the receiver to detect the light scattered back by clouds and precipitation.From the time delay between each outgoing pulse and the backscattered signal, the distance to the scatterer is inferred[43].The VIIRS pixels are, first, identified by the VIIRS Cloud Mask intermediate product, and then the extracted clear-sky pixels are further screened by the field measurements.The pixels with more than one detected cloud base are considered as cloudy, while those with no significant backscatter detected are regarded as true clear-sky pixels.The developed liner models are used to retrieve LWUP.The validation results of LWUP at site C1 are displayed in Figure7.The bias and RMSE of cloudy pixels are −6.65 W/m 2 and 16.04 W/m 2 , while bias and RMSE of clear-sky pixels are −2.94W/m 2 and 15.62 W/m 2 .The VIIRS Cloud Mask intermediate product actually cannot identify all kinds of cloud types, and undetected cloud will cause the underestimation of LWUP.

Figure 7 .
Figure 7. Validation results of the linear model at the Atmospheric Radiation Measurement (ARM) program Southern Great Plains (SGP) C1 site.(a) Pixels that were identified as clear sky by the VIIRS cloud mask whereas cloudy by ground-based Lidar; (b) pixels that were identified as clear sky by both the VIIRS Cloud Mask and the Lidar.

Figure 7 .
Figure 7. Validation results of the linear model at the Atmospheric Radiation Measurement (ARM) program Southern Great Plains (SGP) C1 site.(a) Pixels that were identified as clear sky by the VIIRS cloud mask whereas cloudy by ground-based Lidar; (b) pixels that were identified as clear sky by both the VIIRS Cloud Mask and the Lidar.

Table 1 .
Information of sites in the Surface Radiation Budget Network (SURFRAD) network.

Table 2 .
Fitting results for the linear models and the Multivariate Adaptive Regression Spline (MARS) models.