Radiometric Cross-Calibration of GF-1 PMS Sensor with a New BRDF Model

On-orbit radiometric calibration of a space-borne sensor is of great importance for quantitative remote sensing applications. Cross-calibration is a common method with high calibration accuracy, and the core and emphasis of this method is to select the appropriate reference satellite sensor. As for the cross-calibration of high-spatial resolution and narrow-swath sensor, however, there are some scientific issues, such as large observation angles of reference image, and non-synchronization (or quasi-synchronization) between the imaging date of reference image and the date of sensor to be calibrated, which affects the accuracy of cross-calibration to a certain degree. Therefore, taking the GaoFen-1 (GF-1) Panchromatic and Multi-Spectral (PMS) sensor as an example in this research, an innovative radiometric cross-calibration method is proposed to overcome this bottleneck. Firstly, according a set of criteria, valid MODIS (Moderate Resolution Imagine Spectroradiometer) images of sunny day in one year over the Dunhuang radiometric calibration site in China are extracted, and a new and distinctive bidirectional reflectance distribution function (BRDF) model based on top-of-atmosphere (TOA) reflectance and imaging angles of the sunny day MODIS images is constructed. Subsequently, the cross-calibration of PMS sensor at Dunhuang and Golmud radiation calibration test sites is carried out by using the method presented in this paper, taking the MODIS image with large solar and observation angles and Landsat 8 Operational Land Imager (OLI) with different dates from PMS as reference. The validation results of the calibration coefficients indicate that our proposed method can acquire high calibration accuracy, and the total calibration uncertainties of PMS using MODIS as reference sensor are less than 6%.


Introduction
Radiometric calibration is to convert the digital number (DN) value of a satellite sensor to a physical characterization, such as apparent radiance or apparent reflectance, in the visible and near-infrared sensors.Radiometric calibration is a critical step and the basis of the quantification of remote sensing application.The radiometric calibration not only characterizes the status and degradation with the time of the satellite sensors after launch, but also establishes the physical relationship between multi-source remote sensing data and expands the data source of quantitative remote sensing.Therefore, high precision radiometric calibration plays an important role in the development of quantitative remote sensing.
MODIS bi-directional reflectance distribution function (BRDF) products.Zhong et al. developed a new cross-calibration technique for HJ-1satellite CCD (Charge-coupled Device) sensor [22].In this method, the Landsat ETM+ and Advanced Space-borne Thermal Emission and Reflection Radiometer global digital elevation model (ASTER GDEM) product are used to develop a BRDF model of a desert site.Yang A. et al. updated the method proposed by Zhong et al. and used for the GF-1 satellite WFV sensor [23].The Landsat 8 and the DEM extracted by the three-line camera sensor (TLC) onboard ZiYuan-3 are replaced to develop BRDF model of a desert site.Yang A. et al. also applied this technique to calibrate the GF-4 satellite PMS sensor [24].It can be seen that the BRDF models in prior studies are based on ground measured data or existing BRDF products.Moreover, the MODIS BRDF product and ground reflectance BRDF model are based on the given solar angle, which has some drawbacks and limitations.Therefore, as the absence of ground reflectance and without using BRDF products, how to construct a TOA BRDF model of calibration site and correct the TOA reflectance with a large observation angle (such as MODIS images) to verify the accuracy of calibration results is the main purpose of this paper.
The structure of this paper is as follows: GF-1 PMS, MODIS, Landsat 8 OLI, the test sites, and datasets are described in Section 2. Section 3 introduces the methods of this paper, including BRDF model construction method and radiometric cross-calibration method.The BRDF model construction results, the radiometric cross-calibration results based on BRDF model and the verification of the accuracy of the method proposed in this paper are illustrated in Section 4. Section 5 discusses the uncertainty caused by ground measured data, BRDF model, etc., and the total uncertainty of cross-calibration using MODIS as reference sensor are given.The conclusions of our research are made in Section 6.

Satellites
GF-1 satellite carries six sensors, including four WFV sensors (WFV1, WFV2, WFV3, and WFV4) and two PMS sensors (PMS1 and PMS2).The WFV sensor is a multispectral sensor with four channels, which are blue, green, red, and near infrared (NIR) [25].The spatial resolution of WFV is 16 m and the swath of each WFV sensor is about 200 km, covering 800 km with four WFV sensors.The PMS sensor has four multi-spectral bands with 7.79 m spatial resolution and a panchromatic band with 1.95 m spatial resolution.The swath of each PMS sensor is about 35 km.Both WFV and PMS sensors are linear CCD detectors.The pixel number of multi-spectral and panchromatic band of PMS sensor are 4548 and 18192, respectively.In this paper, just the PMS1 sensor was cross-calibrated using MODIS and OLI.
The MODIS sensor was launched on 18 December 1999 aboard the Terra satellite and has been in operation more than 18 years, providing the global data that are used to monitor long-term changes in the Earth system.The Band 1-4 of MODIS were used to cross-calibrate the PMS1 sensor.The spatial resolutions of MODIS at nadir are 250 m for bands 1 and 2 and 500 m for bands 3 and 4. The reflective solar band (RSB) calibration of MODIS is based on the regular measurements from the solar diffuser (SD) and solar diffuser stability monitor (SDSM), in addition to the monthly scheduled lunar observations.The measurements from SD and SDSM are used for absolute calibration, whereas the lunar measurements along with the SD measurements are used to track the change in on-orbit response vs. scan angle [26,27].In addition to its onboard calibration and characterization activities, MODIS has been calibrated vicariously using Earth targets such as Railroad Valley [28].A 2% uncertainty requirement is specified for the MODIS top-of-atmosphere (TOA) reflectance products, and 5% uncertainty is specified for the MODIS at-sensor spectral radiance products for observations made within ±45 • viewing angles.
Landsat-8, the latest in the Landsat series of satellites, which was launched on 11 February 2013, carries the Operational Land Imager (OLI) as one of its payloads.The OLI is a push broom instrument, whereas all previous recent Landsat instruments are electromechanical (whiskbroom) instruments.The swath of OLI is 185 km and spatial resolution is 30 m for multispectral bands and 15 m for panchromatic band.Two additional multispectral channels were added at 0.445 µm and 1.372 µm for coastal aerosol studies and cirrus cloud detection, respectively.The relative spectral response (RSR) of OLI's NIR band is greatly changed compared with ETM+ to avoid the effect of water absorption in 0.825 µm.MODIS has similar spectral response function compared with OLI.The relative spectral response (RSR) of OLI's panchromatic band is narrower than ETM+ for better distinguishing the vegetation and non-vegetation.The calibration of Landsat 8 OLI is based on inner lamps, solar diffusers, and lunar calibration.The calibration uncertainties of OLI's RSB are about 3% for at-sensor radiance and 2% for at-sensor reflectance [29].
In this paper, the PMS1 sensor is cross-calibrated using both MODIS and OLI.The RSR of these sensors are plotted in Figure 1.The thick black curve is the RSR of PMS, and the red and green thin curves are RSRs of OLI and MODIS, respectively.The RSR of PMS is wider than OLI and MODIS, especially for NIR band.For blue and green bands, the RSR of PMS is more similar with OLI than with MODIS.The RSR of OLI and MODIS in red bands is much closer with each other and narrower than PMS.The center wavelength and solar irradiance of PMS, OLI, and MODIS are summarized in Table 1.The center wavelength is computed by Equation (1).
where τ sensor,i is the RSR of sensor in i-th band (unit: 1).The 'sensor' in this paper represents PMS, OLI, and MODIS.λ is the wavelength (unit: nm).The top of atmosphere (TOA) solar irradiance is computed by Equation ( 2).
where E sensor,i is the TOA solar irradiance of sensor i-th band.E(λ) is the TOA solar irradiance profile range from 400 to 1000 nm (unit: W/m 2 /nm) based on Thuillier solar irradiance model [30].
response (RSR) of OLI's NIR band is greatly changed compared with ETM+ to avoid the effect of water absorption in 0.825 μm.MODIS has similar spectral response function compared with OLI.
The relative spectral response (RSR) of OLI's panchromatic band is narrower than ETM+ for better distinguishing the vegetation and non-vegetation.The calibration of Landsat 8 OLI is based on inner lamps, solar diffusers, and lunar calibration.The calibration uncertainties of OLI's RSB are about 3% for at-sensor radiance and 2% for at-sensor reflectance [29].
In this paper, the PMS1 sensor is cross-calibrated using both MODIS and OLI.The RSR of these sensors are plotted in Figure 1.The thick black curve is the RSR of PMS, and the red and green thin curves are RSRs of OLI and MODIS, respectively.The RSR of PMS is wider than OLI and MODIS, especially for NIR band.For blue and green bands, the RSR of PMS is more similar with OLI than with MODIS.The RSR of OLI and MODIS in red bands is much closer with each other and narrower than PMS.The center wavelength and solar irradiance of PMS, OLI, and MODIS are summarized in Table 1.The center wavelength is computed by Equation (1).
where τ sensor,i is the RSR of sensor in i -th band (unit: 1).The 'sensor' in this paper represents PMS, OLI, and MODIS.λ is the wavelength (unit: nm).
The top of atmosphere (TOA) solar irradiance is computed by Equation (2).
where sensor,i E is the TOA solar irradiance of sensor i-th band.

E
is the TOA solar irradiance profile range from 400 to 1000 nm (unit: W/m 2 /nm) based on Thuillier solar irradiance model [30].E and the altitude is 2894 m.The ground spectrum and atmospheric characteristics of Golmud test site such as aerosol optical thickness and water vapor column are unknown as no ground experiment was carried out at this site before.The homogeneous area of Golmud test site is about 15 km × 15 km. Figure 2 shows two Landsat 8 OLI images over the Dunhuang test site and Golmud test site.The red square boxes in these two images are core test areas with the area of 1.5 km × 1.5 km.An assumption was made that the ground parameters of Dunhuang could be used for Golmud (although at very different elevation and acquired at different times or even days), and the uncertainty due to that assumption will be addressed later in the paper.

Test Sites
The two test sites are used in cross-calibration in this paper, which are the Dunhuang test site and Golmud test site.The Dunhuang test site is an on-orbit radiometric calibration test site for visible/near infrared sensors in China.It has been used for calibrating Chinese satellites since the late 1990s.The Dunhuang test site is large, flat, and homogeneous, located 15 km west of Dunhuang city.It is located at the coordinates 40.1°N and 94.4°E and at an elevation of 1253 m above sea level.The area of the Dunhuang test site is about 30 km × 30 km.There is low probability here of interference from clouds, and low aerosol loading, except during the sandstorm periods in spring.The ground is made up of sand and small black stones with no vegetation.
The Golmud test site is another flat and homogenous test site located 60 km west of Golmud city, Qinghai Province, China.The latitude and longitude of Golmud test site are 36.38°N and 94.23° E and the altitude is 2894 m.The ground spectrum and atmospheric characteristics of Golmud test site such as aerosol optical thickness and water vapor column are unknown as no ground experiment was carried out at this site before.The homogeneous area of Golmud test site is about 15 km × 15 km. Figure 2 shows two Landsat 8 OLI images over the Dunhuang test site and Golmud test site.The red square boxes in these two images are core test areas with the area of 1.5 km × 1.5 km.An assumption was made that the ground parameters of Dunhuang could be used for Golmud (although at very different elevation and acquired at different times or even days), and the uncertainty due to that assumption will be addressed later in the paper.

Datasets
Two image triplets are used for cross-calibration.The first image triplet is collected at the Dunhuang test site.The second image triplet is imaged on Golmud test site.The image acquisition dates and times of PMS, MODIS, and OLI over Dunhuang and Golmud test sites are shown in Table 2.The ground data of Dunhuang test site are used for cross calibration in this paper.The ground reflectance and aerosol optical thickness of the Dunhuang test site were measured on 7 August and 12 August 2014.MODIS and PMS are imaged on 7 August and OLI is imaged on 12 August.The sounding data were also measured at the Dunhuang test site.The spectra vista corporation (SVC) HR1024 spectrometer was used to measure the ground reflectance and CE318 Sun Photometer was used to retrieve the aerosol optical thickness.Figure 3 shows the ground reflectance of the Dunhuang test site, which varied from 0.15 to 0.3.The reflectance is increased with the increase in wavelength.The variance coefficient (the ratio of standard deviation and mean reflectance) of the Dunhuang test site was about 5% as shown in Figure 3.The aerosol optical thickness was 0.367 when PMS1 overpassed Dunhuang test site, with the surface meteorological range of 10.2 km.The water vapor profile of the Dunhuang test site is used in MODTRAN calculation, with the water vapor column of 1.982 gm/cm 2 .For OLI sensor, the surface meteorological range on 12 August was 30.6 km and the water vapor column was 1.243 gm/cm 2 .The ground data of Dunhuang test site are used for cross calibration in this paper.The ground reflectance and aerosol optical thickness of the Dunhuang test site were measured on 7 August and 12 August 2014.MODIS and PMS are imaged on 7 August and OLI is imaged on 12 August.The sounding data were also measured at the Dunhuang test site.The spectra vista corporation (SVC) HR1024 spectrometer was used to measure the ground reflectance and CE318 Sun Photometer was used to retrieve the aerosol optical thickness.Figure 3 shows the ground reflectance of the Dunhuang test site, which varied from 0.15 to 0.3.The reflectance is increased with the increase in wavelength.The variance coefficient (the ratio of standard deviation and mean reflectance) of the Dunhuang test site was about 5% as shown in Figure 3.The aerosol optical thickness was 0.367 when PMS1 overpassed Dunhuang test site, with the surface meteorological range of 10.2 km.The water vapor profile of the Dunhuang test site is used in MODTRAN calculation, with the water vapor column of 1.982 gm/cm 2 .For OLI sensor, the surface meteorological range on 12 August was 30.6 km and the water vapor column was 1.243 gm/cm 2 .As no ground data are available on the Golmud test site, the ground reflectance of the Dunhuang test site is used for SBAF (spectral band adjustment factors) correction.The surface meteorological range of Golmud is assumed to be 10 km and the water vapor column value of 0.2077 gm/cm 2 of mid-latitude winter atmosphere mode is used in MODTRAN calculation.As no ground data are available on the Golmud test site, the ground reflectance of the Dunhuang test site is used for SBAF (spectral band adjustment factors) correction.The surface meteorological range of Golmud is assumed to be 10 km and the water vapor column value of 0.2077 gm/cm 2 of mid-latitude winter atmosphere mode is used in MODTRAN calculation.

Construction Method of BRDF Model
In this paper, a new BRDF correction method based on TOA reflectance of time-series MODIS images is proposed for cross calibration which is mainly used for large-view zenith angle image pairs.Compared with the existing BRDF models, there are several advantages of this new BRDF construction method.Firstly, this method uses TOA reflectance as an input parameter, and the influence of atmospheric orientation as well as the change of surface direction is taken into account.Secondly, the BRDF model in this paper is constructed by using MODIS TOA reflectance data of calibrated test site on sunny days all year round, while MODIS BRDF products are calculated from multi-date, atmospherically corrected data measured over 16-day periods [31].Thirdly, the new BRDF model includes the various solar angles and satellite observation angles from time-series MODIS images of the calibration test sites, rather than based on the given solar angle.The method is based on the following assumptions: For a calibration test site with clear weather and uniform and stable surface, the TOA reflectance at vertical observation angle is assumed to have a high stability throughout the year, and the TOA reflectance differences of different dates are mainly caused by the differences of observation angles and solar angles.Based on these assumptions, a BRDF model of calibration test site based on TOA reflectance is constructed by acquiring MODIS TOA reflectance from different solar and observation angles.One of the difficulties to be solved in this method is how to screen TOA reflectance on sunny days and cloudy days.
A flow chart of the BRDF model construction method is shown in Figure 4, including seven steps.
(1) The time-series MODIS satellite images of the calibration site were collected, and the corresponding DN value, calibration coefficients, and imaging geometry angles were calculated.
where ρ i is the TOA reflectance of i-th MODIS band.The gain re f i and o f f set re f i are the reflectance gain coefficient and reflectance offset coefficient of MODIS, which can be obtained from the MODIS file.The parameter θ s is the solar zenith angle of MODIS image over the test site.The TOA BT of band 31 is calculated by Equation ( 4): where h is Planck constant, k is Boltzmann constant, c is light velocity, λ is wavelength, and L is radiance calculated by Equation ( 5): where, L i is the radiance of i-th band of MODIS.The gain rad i and o f f set rad i are the radiance gain coefficient and radiance offset coefficient of i-th band of MODIS.
(3) The envelope brightness temperature (EBT) is established by TOA BT.In geometry, the envelope of a family of curves is a curve with at least one tangent to each line of the family (the family of curves is the infinite set of some curves, and they have some specific relations.)Then, the image whose difference between EBT and TOA BT is less than 10 K is selected (details can be found in Appendix A).
(4) The variation coefficient (VC) (the ratio of the standard deviation and mean reflectance) of gray value of the calibration site image in the MODIS band 1 was calculated, and the image with VC less than 4% was selected.(5) Eliminate the images with solar zenith angle greater than 55 • .As to the two test sites used in this paper, under the circumstance of the solar zenith angle of MODIS imagery is greater than 55 • , which indicates that the imaging date is in winter, the ground may snow, and the directional reflection will be abnormal.(6) The images which satisfy the difference between EBT and BT less than 10 K, solar zenith angle smaller than 55 • , and VC less than 4% are used as the valid images on a sunny day of the calibration site.Then, the TOA reflectance and imaging angles of the valid calibration site images are calculated.(7) The TOA reflectance BRDF correction coefficients based on the valid calibration site TOA reflectance image and imaging angles are calculated by using the kernel-derived model [32].
The kernel-derived model is represented as follows: where less than 4% was selected.(5) Eliminate the images with solar zenith angle greater than 55°.As to the two test sites used in this paper, under the circumstance of the solar zenith angle of MODIS imagery is greater than 55°, which indicates that the imaging date is in winter, the ground may snow, and the directional reflection will be abnormal.(6) The images which satisfy the difference between EBT and BT less than 10 K, solar zenith angle smaller than 55°, and VC less than 4% are used as the valid images on a sunny day of the calibration site.Then, the TOA reflectance and imaging angles of the valid calibration site images are calculated.(7) The TOA reflectance BRDF correction coefficients based on the valid calibration site TOA reflectance image and imaging angles are calculated by using the kernel-derived model [32].
The kernel-derived model is represented as follows: where θ θ ϕ

BRDF Model Construction
The key of the BRDF model construction method proposed in this paper is to screen the images of sunny day which are extracted by using TOA BT and VC.At the Dunhuang and Golmud sites, a 6

BRDF Model Construction
The key of the BRDF model construction method proposed in this paper is to screen the images of sunny day which are extracted by using TOA BT and VC.At the Dunhuang and Golmud sites, a 6 × 6-pixel (6 km × 6 km) region of interest (ROI) is selected, and the time-series TOA reflectance of the ROI in 2014 is calculated.Then, TOA BT in band 31 and VC in band 1 are calculated by average TOA radiance.It is assumed that the TOA-BT variation of adjacent date images should be less than 10 K in clear weather.Then, we built the EBT and compared the difference between EBT and BT of each image.The image with the temperature difference less than 10 K is regarded as imaged in a sunny day, otherwise, cirrus spissatus existed in the image.Figure 5 shows the extraction results of a sunny day.Moreover, it is assumed that if the VC of the image is more than 4% it is due to the presence of cloud fragmentation.Then, the image is excluded from sunny day images.Figure 6 shows the extraction results of TOA reflectance (MODIS band 1) from sunny day images with VC less than 4% and the difference less than 10 K between TOA BT and EBT.Finally, images satisfying the difference between EBT and TOA BT less than 10 K and VC less than 4% are selected as valid sunny day images, and the information of TOA reflectance and observation angles (including solar zenith angle, solar azimuth angle, view zenith angle, and view azimuth) are extracted.Then, the BRDF correction coefficients are calculated using the kernel-derived model and are shown in Table 3.

Site
Band fiso fgeo fvol Dunhuang Blue 0.2864 0.0525 0.0509 Moreover, it is assumed that if the VC of the image is more than 4% it is due to the presence of cloud fragmentation.Then, the image is excluded from sunny day images.Figure 6 shows the extraction results of TOA reflectance (MODIS band 1) from sunny day images with VC less than 4% and the difference less than 10 K between TOA BT and EBT.Moreover, it is assumed that if the VC of the image is more than 4% it is due to the presence of cloud fragmentation.Then, the image is excluded from sunny day images.Figure 6 shows the extraction results of TOA reflectance (MODIS band 1) from sunny day images with VC less than 4% and the difference less than 10 K between TOA BT and EBT.Finally, images satisfying the difference between EBT and TOA BT less than 10 K and VC less than 4% are selected as valid sunny day images, and the information of TOA reflectance and observation angles (including solar zenith angle, solar azimuth angle, view zenith angle, and view azimuth) are extracted.Then, the BRDF correction coefficients are calculated using the kernel-derived model and are shown in Table 3.

Site
Band fiso fgeo fvol Dunhuang Blue 0.2864 0.0525 0.0509 Finally, images satisfying the difference between EBT and TOA BT less than 10 K and VC less than 4% are selected as valid sunny day images, and the information of TOA reflectance and observation angles (including solar zenith angle, solar azimuth angle, view zenith angle, and view azimuth) are extracted.Then, the BRDF correction coefficients are calculated using the kernel-derived model and are shown in Table 3.

Cross-Calibration Based on BRDF Model
The GF-1 PMS sensor is cross-calibrated in this paper, taking MODIS and OLI as the reference sensor, respectively.
Based on the proposed BRDF model above, the conventional ρ m R TOA reflectance of reference sensor has been replaced by the ρ m,0 R term.
The SBAF (spectral band adjustment factors) was calculated to compensate the spectral response difference, and the SBAF coefficient was defined as Equation ( 8): where ρ m PMS and ρ m R are the measured TOA reflectance for PMS and reference sensor, respectively.These values represent the reflectance directly acquired from multispectral sensor imagery by converting calibrated DNs to TOA reflectance.The ρ s PMS and ρ s R are the simulated TOA reflectance for PMS and reference sensor, respectively.They refer to the reflectance obtained from integrating the RSR of a multispectral sensor with a hyper-spectral profile target.The PMS and reference sensor term can be expressed as Equation ( 9): The hyper-spectral profile ρ(λ) was simulated using a MODTRAN radiative transfer model, according to PMS observation geometry, as well as ground and atmospheric parameters measured over the test site.Then, the SBAF can then be defined as Equation (10): Then, the PMS TOA reflectance over the test site can be calculated as Equation (11): The PMS TOA radiance can be calculated based on solar zenith angle and TOA reflectance according to the Equation ( 12): where L m PMS is the spectral radiance at the sensor's aperture (W/(m 2 sr µm)), d is the Earth-Sun distance measured in astronomical units (AUs), ESUN PMS represents the mean exo-atmospheric solar irradiance (W/(m 2 µm)), and θ s is the solar zenith angle measured in • .
The PMS calibration coefficient gain is defined as Equation ( 13):

Cross-Calibration Results Based on the BRDF Model
According to the proposed BRDF model and cross-calibration method proposed in this paper, the radiometric calibration coefficients of the GF-1 satellite PMS1 sensor with BRDF correction over Dunhuang and Golmud test sites were calculated.Meanwhile, the cross-calibration coefficients without BRDF correction were also calculated for comparison.Table 4 illustrates the BRDF correction coefficients (COE BRDF ), cross-calibration coefficients with and without BRDF correction at Dunhuang and Golmud test sites by using Landsat 8 OLI sensor and MODIS sensor, and the relative error of cross-calibration coefficients between the two sensors.The relative error is defined as "1-cross-calibration coefficient by using Landsat 8 OLI sensor/cross-calibration coefficient by using MODIS sensor".The site calibration coefficients and the relative error between cross-calibration coefficients with and without BRDF correction by using Landsat 8 OLI sensor and MODIS sensor and site calibration coefficients is exhibited in Table 5.The relative error is defined as "1-cross-calibration coefficient/site calibration coefficient". 1'CRESDA' represents the official calibration coefficients by site calibration method. 2 'L' represents the cross-calibration coefficients by using Landsat 8 OLI sensor. 3'M' represents the cross-calibration coefficients by using the MODIS sensor. 4'1-M/SC' represents the relative error between cross-calibration coefficients by using MODIS sensor and site calibration coefficients.

Consistency Analysis of Calibration Results Using Different Test Sites
As no ground data are available on the Golmud test site, the ground reflectance of the Dunhuang test site is used as alternative data source.Therefore, the consistency of the calibration results obtained by using Dunhuang and Golmud sites is analyzed, and the relative errors (1-Dunhuang cross-calibration coefficient/Golmud cross-calibration coefficient) of the calibration results of the two sites are calculated, which are shown in Table 6.It can be seen from the table that the variation of relative error using OLI as reference sensor ranges from 3.22% to 6.10%, and −0.89% to −7.70% using MODIS as a reference sensor.The absolute value of the maximum relative error is 7.70% over the green band by using MODIS as reference sensor.This indicates that the cross-calibration results of the two test sites have good consistency.It is proved that the Golmud site can be used for radiometric calibration of remote sensing satellites even without ground-based measured data.

Consistency Analysis of Calibration Results Using Different Reference Satellites
Table 4 shows the results of cross-calibration of GF-1 satellite PMS1 sensor using MODIS and Landsat 8 OLI as reference sensors, respectively.The relative errors (1-cross-calibration coefficient using Landsat 8 OLI as a reference/cross-calibration coefficient using MODIS as a reference) of cross-calibration coefficients between the two reference sensors are calculated.As can be seen from Table 4, the absolute value of relative error of the Golmud calibration test site is less than that of the Dunhuang calibration test site.At the Dunhuang test site, the relative errors of cross-calibration coefficients between the two reference sensors are 5.80%, 8.53%, 6.20%, and 5.87%, while at the Golmud site, the relative errors are −1.08%,−2.66%, −3.17%, and −1.47%, respectively, which means that whether at the Dunhuang or Golmud test site, the cross-calibration results using the two reference satellites are in good agreement, especially at the Golmud site.The results demonstrate that the BRDF-corrected large observation angle MODIS images can be used as reference data for remote sensing radiometric cross-calibration.

Validation of Cross-Calibration Coefficients Accuracy
In order to validate the radiometric cross-calibration coefficients accuracy of the proposed method in this paper, the calibration coefficients with BRDF correction, calibration coefficients without BRDF correction, and the site calibration coefficients are compared, and the relative errors between them are calculated.The relative errors are defined as "1-cross-calibration coefficients/site calibration coefficients" and listed in Table 5.
It can be seen that, obviously, the relative error decreases significantly after BRDF correction.Particularly, when MODIS sensor images are used as reference, the relative errors of each band become lower from 14.06%, 10.99%, 10.85%, and 4.59% without BRDF correction to −2.47%, −4.43%, −2.26%, and −4.76% with BRDF correction over the Dunhuang test site, while at the Golmud test site, from −12.02%, −7.29%, −0.69%, and −7.84% without BRDF correction to 0.36%, 3.04%, 3.64%, and −3.83% with BRDF correction.This is because the MODIS images of calibration test sites have large sensor zenith angles-49.68• of the Dunhuang test site and 53.12 • of the Golmud test site.
On the whole, the absolute value of relative errors of cross-calibration results obtained at both Dunhuang and Golmud test sites are all smaller than 5%, whether Landsat 8 OLI or BRDF-corrected MODIS images are used as a reference.Furthermore, after BRDF correction, the relative error of the Golmud test site is smaller than that of the Dunhuang test site for most bands, except for the NIR band (using Landsat 8 OLI as reference sensor) and the red band (using MODIS as reference sensor).In general, the relative errors illustrate that the cross-calibration method proposed in this paper has high radiometric calibration accuracy.

Discussion
Normalized MODIS TOA reflectance takes the place of the original MODIS TOA reflectance in the proposed method.Taking the cross-calibration of PMS using MODIS reference sensor over the Golmud test site as an example, the uncertainty of this paper is analyzed as below.

Uncertainty Caused by Ground-Measured Parameters
Ground-measured parameters include ground spectra and atmospheric parameters.As no ground-measured spectral and atmospheric parameters of the Golmud test site are available in this paper, it was assumed that the Golmud test site has the same spectral and atmospheric characteristics as Dunhuang test site, which will affect the accuracy of SBAF.Consequently, in order to analyze the influence of ground measured parameters, new sand spectrum, atmospheric parameters, and aerosol models were used to calculate new SBAF and new cross-calibration coefficients from the Golmud test site.The Dunhuang ground-measured spectrum was replaced by the typical sand spectrum (as shown in Figure 7).The aerosol type was set to be rural type, the surface meteorological range of PMS was modified from 10 KM to 23 KM (a default value in MODTRAN software when the aerosol type is set to rural type) and the water vapor column was changed from 0.2077 gm/cm 2 to 3.00 gm/cm 2 .Then, the new cross-calibration coefficients were calculated, and the uncertainty was analyzed, as shown in Table 7.The results show that, among the four influence factors, the site ground spectrum had the greatest impact, with the first band at 4.79%, and the other three bands at 0.67%, 1.05%, and −1.40%, respectively.This is because the RSRs of PMS and MOIS in blue bands have larger difference than other bands.The total uncertainty of SBAF is 4.85% in the first band and less than 2% in the other three bands.

Uncertainty Caused by BRDF Model
The BRDF model proposed in this paper is based on time series of MODIS TOA reflectance.In order to analyze the accuracy of BRDF model, a MODIS image of Golmud test site on 24 February 2014 is selected, and DN value and imaging geometry information of the site are extracted.Then, the official reflectance calibration coefficients of MODIS and BRDF model coefficients are used to calculate the measured TOA reflectance and simulated TOA reflectance, respectively, and the relative error between them is calculated.The results are summarized in Table 8.As can be seen from the table, the relative errors of the four bands are −0.60%,−1.27%, −1.16% and 0.97%, respectively, which illustrates that the BRDF model constructed in this paper has high accuracy.

Total Uncertainty of Cross-Calibration Using MODIS as Reference Sensor
The total uncertainty of cross calibration is summarized in this part.There are five influence factors in cross-calibration.
(1) SBAF: In some cases, ground-measured parameters such as spectral reflectance, surface meteorological range, and water vapor column are not available for most cross calibration.Therefore, empirical ground-measured parameters could be used in cross calibration.In this article, the ground measured parameters of Dunhuang are used to calculate the SBAF coefficient of Golmud.The uncertainty caused by SBAF was 4.85%, 0.67%, 1.06%, and 1.59%, respectively.

Uncertainty Caused by BRDF Model
The BRDF model proposed in this paper is based on time series of MODIS TOA reflectance.In order to analyze the accuracy of BRDF model, a MODIS image of Golmud test site on 24 February 2014 is selected, and DN value and imaging geometry information of the site are extracted.Then, the official reflectance calibration coefficients of MODIS and BRDF model coefficients are used to calculate the measured TOA reflectance and simulated TOA reflectance, respectively, and the relative error between them is calculated.The results are summarized in Table 8.As can be seen from the table, the relative errors of the four bands are −0.60%,−1.27%, −1.16% and 0.97%, respectively, which illustrates that the BRDF model constructed in this paper has high accuracy.

Total Uncertainty of Cross-Calibration Using MODIS as Reference Sensor
The total uncertainty of cross calibration is summarized in this part.There are five influence factors in cross-calibration.
(1) SBAF: In some cases, ground-measured parameters such as spectral reflectance, surface meteorological range, and water vapor column are not available for most cross calibration.Therefore, empirical ground-measured parameters could be used in cross calibration.In this article, the ground measured parameters of Dunhuang are used to calculate the SBAF coefficient of Golmud.The uncertainty caused by SBAF was 4.85%, 0.67%, 1.06%, and 1.59%, respectively.
(2) BRDF model: The construction of BRDF model in this paper is based on time-series of MODIS image data and has high accuracy.The uncertainty of BRDF model is −0.60%, −1.27%, −1.16%, and 0.97%, respectively.(3) Radiative transfer model: The simulation software MODTRAN is used in this paper, which might have some effects on the TOA reflectance, and the uncertainty of MODTRAN is about 1%.(4) Reference sensor: MODIS and Landsat 8 OLI were used as reference sensor in this paper.The uncertainty of the MODIS and OLI on board calibration was within 2%, which was the basic uncertainty for the PMS1 cross-calibration.( 5) Image uncertainty: The image uncertainty mainly comes from geometric mismatching between PMS image and reference image.In our study, the DNs were mean value of homogenous areas, so the image uncertainty of geometric mismatching was within 0.5%.
Table 9 shows the uncertainties caused by various factors and the total uncertainties.The total uncertainties are calculated by square root methods, which are 5.40%, 2.70%, 2.78%, and 2.95% of four bands, respectively.

Results and Discussion of Time Series Cross-Calibration
For further study, the PMS1 data on the Dunhuang and Golmud test sites since the launch of GF-1 satellite were collected, and we also obtained the corresponding MODIS and Landsat 8 OLI data.However, as the PMS1 gain was adjusted by the CRESDA after September 2016, the subsequent data are not comparable.Consequently, from the launch date of GF-1 satellite to September 2016, there are another five valid PMS1 data that has been used in research.Among them, there are two scenes at Dunhuang test site and three scenes on Golmud test site.The cross-calibration of PMS1 on the two test sites (Dunhuang, Golmud) is carried out by using the method presented in this paper.The calibration results were shown in Table 10.The relative error between cross-calibration coefficients and CRESDA site calibration coefficients were shown in Table 11.The relative errors are defined as "1-cross-calibration coefficients/site calibration coefficients".It can be seen from Tables 10 and 11 that the calibration coefficients obtained by this method were in good agreement, whether at Dunhuang or Golmud sites.Obviously, except for individual bands (NIR band at Golmud test site using MODIS as reference sensor on 12 December 2014), all relative differences were less than 10%, and most of the relative errors were less than 5%.
The relative errors of calibration coefficients between the two reference sensors (Landsat 8 OLI and MODIS) were also calculated, as shown in Table 12.The relative error is defined as "1-cross-calibration coefficient by using Landsat 8 OLI sensor/cross-calibration coefficient by using MODIS sensor".As can be seen from Table 12, the absolute value of relative error of Golmud calibration test site was less than that of Dunhuang calibration test site.In general, the results obtained from different reference sensors were of obvious consistence.The relative errors at the Golmud test site were less than 5%.While, at the Dunhuang test site, most relative errors were less than 9%.In addition, in order to further study the variation of cross-calibration coefficients with time, the variation characteristics of calibration coefficients in the blue band (using MODIS as reference sensor) are drawn, as shown in Figure 8.As can be seen, the GF-1 PMS1 sensor has a certain attenuation, but it is not particularly obvious.The attenuation is −3.39 × 10 −6 per day.However, due to the small amount of valid data to fully explain and discuss the attenuation trend and characteristics of GF-1 PMS1, adequate data needs to be collected from more test sites to complete the issue.
The average values of cross-calibration coefficients using two reference sensors were calculated.Moreover, the relative errors between the average cross-calibration coefficients and the laboratory calibration coefficients and the official site calibration coefficients were also compared.
attenuation, but it is not particularly obvious.The attenuation is −3.39 × 10 −6 per day.However, due to the small amount of valid data to fully explain and discuss the attenuation trend and characteristics of GF-1 PMS1, adequate data needs to be collected from more test sites to complete the issue.The average values of cross-calibration coefficients using two reference sensors were calculated.Moreover, the relative errors between the average cross-calibration coefficients and the laboratory calibration coefficients and the official site calibration coefficients were also compared.
As shown in Table 13, the absolute values of relative errors between average cross-calibration coefficients and laboratory calibration coefficients were between 5% and 18%, which means that PMS1 was attenuated obviously after satellite launching, and the laboratory calibration coefficients was no longer available and post-launch calibration of PMS1 sensor was necessary.However, the absolute values of relative errors between average cross-calibration coefficients and CRESDA site calibration coefficients were less than 3%, except for individual bands. 1 'Lab' represents the calibration coefficients obtained in laboratory conditions. 2 'CRESDA' represents the official calibration coefficients by site calibration method. 3'L' represents the average cross-calibration coefficients by using Landsat 8 OLI as reference sensor. 4'M' represents the average cross-calibration coefficients by using MODIS as reference sensor. 5'(1-L/Lab)%' represents the relative error between average cross-calibration coefficients by using Landsat 8 OLI as reference sensor and laboratory calibration coefficients.
In summary, the results of this method are in good agreement with those of site calibration and can be used to validate the results of site calibration, with higher calibration frequency and lower As shown in Table 13, the absolute values of relative errors between average cross-calibration coefficients and laboratory calibration coefficients were between 5% and 18%, which means that PMS1 was attenuated obviously after satellite launching, and the laboratory calibration coefficients was no longer available and post-launch calibration of PMS1 sensor was necessary.However, the absolute values of relative errors between average cross-calibration coefficients and CRESDA site calibration coefficients were less than 3%, except for individual bands. 1 'Lab' represents the calibration coefficients obtained in laboratory conditions. 2 'CRESDA' represents the official calibration coefficients by site calibration method. 3'L' represents the average cross-calibration coefficients by using Landsat 8 OLI as reference sensor. 4'M' represents the average cross-calibration coefficients by using MODIS as reference sensor. 5'(1-L/Lab)%' represents the relative error between average cross-calibration coefficients by using Landsat 8 OLI as reference sensor and laboratory calibration coefficients.
In summary, the results of this method are in good agreement with those of site calibration and can be used to validate the results of site calibration, with higher calibration frequency and lower cost.

Conclusions
An innovative cross-radiation calibration method is proposed in this paper, which is used to cross-calibrate GF-1 PMS camera at Dunhuang and Golmud test site, respectively, taking MODIS and Landsat 8 OLI as reference sensors.
Firstly, a BRDF model is constructed based on the TOA reflectance of MODIS in long time-series sunny day, which are extracted by the TOA BT and the variation coefficient of TOA reflectance.This process is an important step of this paper.Then the TOA reflectance of MODIS images with large imaging angles is corrected by BRDF model.After that, the radiometric cross-calibration of PMS is carried out by using BRDF-corrected MODIS and OLI images, and the radiometric calibration coefficients are calculated by the least square method.
Compared with previous studies, it is found that the accuracy of cross-radiation calibration has been significantly improved after directional correction of large observation angle MODIS images using the BRDF model constructed in this paper.Furthermore, the study also finds that although ground measured parameters of Dunhuang test site are used as substitution for the Golmud test site in cross-calibration, the calibration coefficients obtained are in good agreement with that of the Dunhuang site.By comparing with the site calibration coefficients and analyzing the uncertainty of this paper, the total calibration uncertainties of PMS using MODIS as reference sensor are about 2%-6%, which are similar with the accuracy of vicarious calibration.The is promising, and in the future, the robustness and universality of the technique proposed in this paper should also be assessed by applying it to other radiometric calibration test sites and similar optical satellite sensors.In addition, as CRESDA has adjusted the gain of GF-1 PMS1 on 25 September 2016, this paper only studies the PMS1 data from the launch of GF-1 satellite on 26 April 2013 to 25 September 2016.The radiation performance and calibration accuracy of the gain-adjusted GF-1 PMS1 sensor will be studied in the follow-up study.(2) A scatter diagram is drawn by DOY as X-axis and BT as Y-axis.As shown in Figure A1.(3) Computing Principle of EBT Firstly, take the first point as the starting point, connect it with the first adjacent point on its right side, and draw an extension line.As shown in Figure A2a.Secondly, determine whether other points fall on the right side of the extension line.If a point falls to the left side of the extension line, the adjacent point is discarded.Then, connect the starting point with the next adjacent point, draw the extension line, and determine again whether all other points fall on the right side of the extension line, and so on, until all points fall on the right side of the extension line.
At this point, the adjacent point is used as a new starting point, and the first adjacent point on its right side is connected, and the extension line is drawn to determine whether the other remaining points are on the right side of the extension line again.
Thirdly, the envelope of bright temperature (EBT) can be obtained by repeating the above process, as shown in Figure A2a.
Finally, the points on the EBT are represented by red solid boxes, as shown in Figure A2b.Secondly, the difference ∆t (∆t = BTD − BTA) between the ordinates (BT) of point A and point D is calculated.If ∆t is less than 10 K, then point A is considered to be sunny day, and represented by a red solid box.Otherwise, if ∆t is greater than 10 K, then this point is not sunny day.
Finally, according to the above principles, other points within the envelope are judged separately.Then, the sunny day that meet the requirement (∆t < 10 K) can be screened out.Secondly, determine whether other points fall on the right side of the extension line.If a point falls to the left side of the extension line, the adjacent point is discarded.Then, connect the starting point with the next adjacent point, draw the extension line, and determine again whether all other points fall on the right side of the extension line, and so on, until all points fall on the right side of the extension line.
At this point, the adjacent point is used as a new starting point, and the first adjacent point on its right side is connected, and the extension line is drawn to determine whether the other remaining points are on the right side of the extension line again.
Thirdly, the envelope of bright temperature (EBT) can be obtained by repeating the above process, as shown in Figure A2a.
Finally, the points on the EBT are represented by red solid boxes, as shown in Figure A2b.Secondly, the difference ∆t (∆t = BTD − BTA) between the ordinates (BT) of point A and point D is calculated.If ∆t is less than 10 K, then point A is considered to be sunny day, and represented by a red solid box.Otherwise, if ∆t is greater than 10 K, then this point is not sunny day.
Finally, according to the above principles, other points within the envelope are judged separately.Then, the sunny day that meet the requirement (∆t < 10 K) can be screened out.
. The Dunhuang test site is an on-orbit radiometric calibration test site for visible/near infrared sensors in China.It has been used for calibrating Chinese satellites since the late 1990s.The Dunhuang test site is large, flat, and homogeneous, located 15 km west of Dunhuang city.It is located at the coordinates 40.1 • N and 94.4 • E and at an elevation of 1253 m above sea level.The area of the Dunhuang test site is about 30 km × 30 km.There is low probability here of interference from clouds, and low aerosol loading, except during the sandstorm periods in spring.The ground is made up of sand and small black stones with no vegetation.The Golmud test site is another flat and homogenous test site located 60 km west of Golmud city, Qinghai Province, China.The latitude and longitude of Golmud test site are 36.38• N and 94.

Figure 2 .
Figure 2. Landsat 8 OLI image over two test sites.The image of Dunhuang test site is imaged on 12 August 2014, and the image of Golmud test site is acquired on 26 February 2014.(a) Dunhuang test site; (b) Golmud test site

Figure 2 .
Figure 2. Landsat 8 OLI image over two test sites.The image of Dunhuang test site is imaged on 12 August 2014, and the image of Golmud test site is acquired on 26 February 2014.(a) Dunhuang test site; (b) Golmud test site.

Figure 3 .
Figure 3. Ground reflectance and variance coefficient of the Dunhuang test site.

Figure 3 .
Figure 3. Ground reflectance and variance coefficient of the Dunhuang test site.

( 2 )
The time series of MODIS TOA reflectance of bands 1-4 and TOA bright temperature (BT) of band 31 are calculated based on calibration coefficients and BT calculation equation.The TOA reflectance of bands 1-4 are calculated by Equation (3):

Figure 5 .
Figure 5. Top-of-atmosphere (TOA) bright temperature (BT) of all images (blue hollow rhombus) and images of sunny days with difference less than 10 K between TOA BT and envelope brightness temperature (EBT) (red solid frame): (a) Dunhuang; (b) Golmud.

Figure 5 .
Figure 5. Top-of-atmosphere (TOA) bright temperature (BT) of all images (blue hollow rhombus) and images of sunny days with difference less than 10 K between TOA BT and envelope brightness temperature (EBT) (red solid frame): (a) Dunhuang; (b) Golmud.

Figure 7 .
Figure 7. Spectra used for SBAF calculation, where Dunhuang is raw spectrum and sand is new spectrum.

Figure 8 .
Figure 8. Time-dependent variation of cross-calibration coefficients in GF-1 PMS1 Blue band.'Lab' represents the laboratory calibration coefficients of blue band; 'CRESDA' represents the site calibration coefficients of blue band.

Figure 8 .
Figure 8. Time-dependent variation of cross-calibration coefficients in GF-1 PMS1 Blue band.'Lab' represents the laboratory calibration coefficients of blue band; 'CRESDA' represents the site calibration coefficients of blue band.

Remote 22 ( 2 )
Sens. 2019, 11,  x FOR PEER REVIEW 19 of A scatter diagram is drawn by DOY as X-axis and BT as Y-axis.As shown in FigureA1.

22 ( 2 )
Remote Sens. 2019, 11, x FOR PEER REVIEW 19 of A scatter diagram is drawn by DOY as X-axis and BT as Y-axis.As shown in Figure A1.

Figure A2 .
Figure A2.Calculation process of envelope of bright temperature (EBT).(a) Process of EBT calculation; (b) Results of EBT calculation.

( 4 )
Principle of sunny day data screeningAs shown in FigureA3, taking point A as an example.Firstly, the perpendicular line is made through point A, intersects with envelope BC at point D, and calculates the ordinate value of point D (which can be calculated by the coordinates of the three points of A, B, and C).

22 Figure A2 .
Figure A2.Calculation process of envelope of bright temperature (EBT).(a)Process of EBT calculation; (b)Results of EBT calculation.

( 4 )
Principle of sunny day data screeningAs shown in FigureA3, taking point A as an example.Firstly, the perpendicular line is made through point A, intersects with envelope BC at point D, and calculates the ordinate value of point D (which can be calculated by the coordinates of the three points of A, B, and C).

Figure A3 .
Figure A3.Map of sunny day screening.

Figure A3 .
Figure A3.Map of sunny day screening.

Table 1 .
Band number, center wavelength, spectral range, and solar irradiance of PMS, MODIS, and OLI.

Table 1 .
Band number, center wavelength, spectral range, and solar irradiance of PMS, MODIS, and OLI.
2.2.Test SitesThe two test sites are used in cross-calibration in this paper, which are the Dunhuang test site and Golmud test site

Table 2 .
Input parameters of MODTRAN for spectral band adjustment factors (SBAF) calculation.

Table 2 .
Input parameters of MODTRAN for spectral band adjustment factors (SBAF) calculation.

Table 3 .
BRDF model parameters based on the MODIS TOA reflectance over the Dunhuang and Golmud test site.

Table 3 .
BRDF model parameters based on the MODIS TOA reflectance over the Dunhuang and Golmud test site.

Table 3 .
BRDF model parameters based on the MODIS TOA reflectance over the Dunhuang and Golmud test site.

Table 4 .
BRDF correction coefficients, radiometric cross-calibration coefficients with and without BRDF correction at Dunhuang and Golmud test sites by using Landsat 8 OLI sensor and MODIS sensor, and the relative error of cross-calibration coefficients between the two sensors.

Table 5 .
Site calibration coefficients and the relative error between cross-calibration coefficients and site calibration coefficients.

Table 6 .
The relative error of cross-calibration coefficients between Dunhuang and Golmud sites using Landsat 8 OLI and MODIS sensor as a reference.
1 'D' represents the Dunhuang site.2'G'represents the Golmud site.3'1-D/G'represents the relative error of cross-calibration coefficients between Dunhuang and Golmud site.

Table 7 .
Remote Sens. 2019, 11,x FOR PEER REVIEW 14 of 22 Spectra used for SBAF calculation, where Dunhuang is raw spectrum and sand is new spectrum.Total uncertainty caused by ground spectrum, atmosphere, and aerosol type.

Table 8 .
Relative error between simulated TOA reflectance calculated by BRDF and measured TOA reflectance calculated by official calibration coefficients.

Table 7 .
Total uncertainty caused by ground spectrum, atmosphere, and aerosol type.

Table 8 .
Relative error between simulated TOA reflectance calculated by BRDF and measured TOA reflectance calculated by official calibration coefficients.

Table 9 .
Total calibration uncertainty of cross-calibration using MODSI over Golmud test site.

Table 10 .
Laboratory calibration coefficients, CRESDA site calibration coefficients and cross-calibration coefficients calculated in this paper.

Table 11 .
The relative errors between cross-calibration coefficients and CRESDA site calibration coefficients.

Table 12 .
The relative errors of cross-calibration coefficients between the two reference sensors.

Table 13 .
Laboratory calibration coefficients, site calibration coefficients and average cross-calibration coefficients and relative error between them.

Table 13 .
Laboratory calibration coefficients, site calibration coefficients and average cross-calibration coefficients and relative error between them.

Table A1 .
Day of year (DOY) and Bright Temperature (BT).