Using HJ-1 CCD and MODIS Fusion Data to Invert HJ-1 NBAR for Time Series Analysis, a Case Study in the Mountain Valley of North China

: HJ-1 charge-coupled device (CCD) data with high temporal and medium spatial resolution are widely used in environmental and disaster monitoring in China. However, due to bad weather, it is difﬁcult to obtain sufﬁcient time-continuous HJ-1 CCD data for environmental monitoring. In this study, the mountain valley with farmland and forestland in North China is selected as the experimental area, and HJ-1 CCD and moderate resolution imaging spectroradiometer (MODIS) data are used in the case study. An improved method of fusing data and inversing surface reﬂectivity is presented to obtain the HJ-1 inversion network-based application resolution (NBAR) data using linear matching of the Ross Thick-Li Sparse Reciprocal (RTLSR) model, and then predicted reﬂectivity using the seasonal autoregressive integrated moving average (SARIMA) model. The fusion data have advantages of high spatial and temporal resolution, as well as meeting the requirements of high quality and quantity of small-scale regional data. This case study provides a feasibility method for the HJ-1 satellites to produce the secondary products for small-scale remote sensing ground surface research. It also provides a reference for dynamic information acquisition and application of small satellite data, contributing to the improvement in RS estimation of surface environment variables.


Introduction
HJ-1 satellites, with high spatial and moderate temporal and spectral resolution data, are the first small satellite constellations dedicated to environmental and disaster monitoring and forecasting in China [1,2].HJ-1 charge-coupled device (CCD) data are widely used in land use classification [3], and also used in wetland classification based on objectoriented classification methods, resulting in a low-cost and effective technical means to obtain high-precision wetland distribution information [4].In other applications, the improved Carnegie-Ames-Stanford Approach (CASA) model is used to estimate the net primary productivity of the vegetation [5,6], and the soil organic matter content in the experimental area was previously monitored by HJ-1 data [7][8][9].
However, due to bad weather, the limited quantity restricts the quality of remote sensing (RS) extraction, the quality of HJ-1 two-level products is poor, geometric calibration is not accurate, and available data are reduced.HJ-1 CCD data are often used in combination with other RS data [10].Fusing with multi-sensor observations is an effective way to improve the observation frequency [11,12].The leaf area index (LAI) is usually calculated using the multi-sensor observation data set, which is constructed by sensors of HJ-1 CCD and Landsat 8 OLI data, and results have shown that the fusion data set can produce LAI products with reliable precision and continuous time resolution [13,14].
The temporal characteristics of surface reflectivity is a key factor for land surface process research, and has a great significance for the prediction of future surface reflectivity, Appl.Sci.2022, 12, 12233 2 of 18 such as hydrological forecasting, crop pest and disease disaster prediction, environmental pollution control, and ecological balance [15].Compared to traditional surface reflectivity, the nadir bidirectional reflectance distribution function (BRDF)-adjusted reflectivity (NBAR) has advantages of being able to provide comparable dynamic data, thus avoiding the angular impact, and reflecting the dynamic variation patterns of the vegetated target area [16,17].Li et al. provided an effective method to simulate the temporal changes in the MODIS NBAR time series of typical farmland surfaces by means of the season-trend statistics [18], and they also used the improved NBAR data as matching values for the LAI neural network estimation, making the predicted LAI time series more continuous than the MODIS LAI [14,19].
In this study, HJ-1 CCD and MODIS data were selected as data sources to invert the HJ-1 NBAR with 30 m resolution based on the MODIS NBAR production algorithm.The inversion data have advantages of high spatial and temporal resolution, as well as meeting the requirements of high quality and quantity of small-scale regional data.This provides a feasibility method for the HJ-1 satellites to produce the secondary products for small-scale RS ground surface research, and provides a reference for dynamic information acquisition and application of small satellite data, contributing to the improvement in RS estimation of surface environment variables.Moreover, the NBAR time series is a special sequence of surface observations, which has a pronounced seasonal change.It contains an important information about the four seasons and solar radiation; its predictions can contribute to early warning of pests and phonological changes, and it can be useful in other studies.

Study Areas
The study area is located in Huailai County in northern Hebei Province of North China, adjacent to Beijing in the East (See Figure 1a for location).The central position is about 115 • east longitude and 40 • north latitude, with a total area of 1801 km 2 .Huailai County is located in the north of the Yanshan Mountains and the upper reaches of Yongding River.It is located in the semi-arid area of the middle temperate zone and belongs to the temperate continental monsoon climate.

Materials
The satellite data used in this study are listed in Table 1: (1) HJ-1A/HJ-1B data: the two-level system geometric correction products of HJ-1A/HJ-1B satellites, which have spatial resolution of 30 m, a revisit period of 4 days, viewing swath width of 700 km, and include a total of four bands (blue, green, red, and nearinfrared).( 2) MODIS NBAR products: the MODIS NBAR products (MOD43A4), which have spatial resolution of 500 m, three-level terrestrial standard, and 16 days of synthetic data Two interesting regions (2 × 2 km) with eight measuring sites were selected as the research areas (see Figure 1b for locations).The north region (40.

Materials
The satellite data used in this study are listed in Table 1: (1) HJ-1A/HJ-1B data: the two-level system geometric correction products of HJ-1A/HJ-1B satellites, which have spatial resolution of 30 m, a revisit period of 4 days, viewing swath width of 700 km, and include a total of four bands (blue, green, red, and near-infrared).( 2) MODIS NBAR products: the MODIS NBAR products (MOD43A4), which have spatial resolution of 500 m, three-level terrestrial standard, and 16 days of synthetic data products.
(3) MODIS daily products: the MODIS global daily reflectivity products (MOD09GA), for which the spatial resolution is 500 m.(4) The surface measured data: the measured reflectance data collected at the same time as the test sites, measured by ASD spectrometer.The red (R) and near-infrared (NIR) bands of the NBAR data were chosen for the time series analysis, as they were sensitive to the vegetation growth and were typically used in most quantifications of vegetation inversion, such as LAI and NDVI, among others.

HJ-1 Data Preprocessing
To obtain more accurate data, we used the earth resources data analysis system (ERDAS) software, based on the Landsat 8 data of the same region, to georeference the HJ-1 data again.Based on the size of the target object and filed survey, 20 control points for geometric correction were selected through a visual recognition method.By selecting the ground object points with obvious visual features on the image, such as road intersections and waterway junction lines, we ensured that the control points were evenly distributed on the two images.
The results showed that the geometric correction only caused shift correction, and there were few changes before and after image correction.The accuracy of Chinese satellite data was gradually increased.
Due to the influence of the atmosphere, product data with geometric corrections were distorted when analyzing the true surface reflectivity.To obtain the real surface information from satellite images, the next step was to accurately remove the atmospheric effect, or atmospheric correction, which eliminates the influence of atmospheric aerosol on clutter reflections and obtains the real surface reflectivity.
To improve the efficiency of atmospheric correction, we calculated the lookup table offline based on the 6S model [20].The table contained different atmospheric aerosol optical thickness, solar zenith angle, observed zenith angle, and relative azimuth angle.It provided the atmospheric correction coefficient to achieve atmospheric correction of HJ-1 A/B data.To omit steps, we used the absolute atmospheric correction methods, while removing the effect of atmosphere in the image, to convert digital number (DN) value into surface reflectivity.
In the 6S model, the atmospheric correction is based on the following equation: where x a , x b , and x c represent the three atmospheric correction coefficients; L i represents the radiance value measured for the i-th band of the sensor; ρ represents the surface reflectance corrected by atmosphere.
When entering a parameter, we selected the predefined standard mode in 6S as the input.The two selected atmospheric models were mid-latitude summer and mid-latitude winter.The aerosol mode selected was the continental type, and assumed that the surface has uniform Lambertian reflection characteristics.
The aerosol optical thickness at 550 nm was set to 0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.8, 1.0, 1.5, and 2.0.Six sensor zeniths were set up within the viewing angle of the CCD camera.Nineteen relative azimuth angles were set up within the range of 0-85 • .Eighteen solar zenith angles were set up within the range of 0-180 • [21].After atmospheric correction, we obtained more realistic surface information, as indicated in Figure 2a.
After removing the atmospheric effect, radiometric calibration of the data was carried out to convert the DN value into the surface reflectivity, based on the spectral response function of the HJ and MODIS data.The four-band spectral response functions of HJ-1A CCD1 and MODIS are shown in Figure 2b.
offline based on the 6S model [20].The table contained different atmospheric aerosol optical thickness, solar zenith angle, observed zenith angle, and relative azimuth angle.It provided the atmospheric correction coefficient to achieve atmospheric correction of HJ-1 A/B data.To omit steps, we used the absolute atmospheric correction methods, while removing the effect of atmosphere in the image, to convert digital number (DN) value into surface reflectivity.
In the 6S model, the atmospheric correction is based on the following equation: where  ,  , and  represent the three atmospheric correction coefficients;  represents the radiance value measured for the i-th band of the sensor;  represents the surface reflectance corrected by atmosphere.
When entering a parameter, we selected the predefined standard mode in 6S as the input.The two selected atmospheric models were mid-latitude summer and mid-latitude winter.The aerosol mode selected was the continental type, and assumed that the surface has uniform Lambertian reflection characteristics.
The aerosol optical thickness at 550 nm was set to 0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.8, 1.0, 1.5, and 2.0.Six sensor zeniths were set up within the viewing angle of the CCD camera.Nineteen relative azimuth angles were set up within the range of 0-85°.Eighteen solar zenith angles were set up within the range of 0-180° [21].After atmospheric correction, we obtained more realistic surface information, as indicated in Figure 2a.
After removing the atmospheric effect, radiometric calibration of the data was carried out to convert the DN value into the surface reflectivity, based on the spectral response function of the HJ and MODIS data.The four-band spectral response functions of HJ-1A CCD1 and MODIS are shown in Figure 2b.

MODIS Data Preprocessing
The coordinate system of MODIS data is different from the conventional ellipsoid system, which uses a regular grid as the imaging unit and has deformation in the longitude and latitude diagonal direction.Therefore, the projection of MODIS data should be transformed into the HJ-1 projection coordinate system of UTM_zone_50N, and then the second geometric correction adjustment on the MODIS data is conducted based on Landsat 8 data, causing the MODIS data and HJ-1 data to overlap accurately.

Measured Data Preprocessing
To measure the surface reflectivity, we used the ASD spectrometer.The ASD spectrometer has a wavelength range of 350-2500 nm, and has reference reflectivity plates of 40%, 50%, and 99%.We set the ASD spectrometer to vertical.To eliminate the effects of random noise, we used the mean surface measured data for 20 min before and after the HJ-1 transit, and collected the solar zenith angles at the same time.
Through the RTLSR model, we calculated the NBAR values, and the average NBAR values at these eight sample points in two interest areas, to evaluate the overall performance.

Method
An overview of the methods is shown in Figure 3, and each step is described in detail in this section.
We focused on inversing HJ-1 NBAR, based on the MOD43 product inversion method, and analysis of time series characteristics of two types of vegetation land cover: forest and farmland.The change in the vegetation within such areas can be well depicted by the time series using HJ-1 NBAR data, which has less mixed pixels than MODIS.
The red and NIR bands of HJ-1 reflectivity data and MOD09GA data were used as input data in the process.We created the HJ-1 fusion data of 16 days as a group, and 8 days as a cycle.In each cycle, it was assumed that the change in the object matter was negligible and the HJ-1 NBAR was reversed.
Using the inversion of the NBAR and the MODIS NBAR data, we conducted a time series feature analysis in both the forest and the farmland areas of interest.Through the prediction and comparison with surface measured values, we proved the practicality of the data.
However, referring to the MODIS flags file, which on behalf of the image point were cloud, cloud shadow, ice, and water, and 15 other types of coverage, we excluded unavailable data from HJ-1 data, resulting in insufficient data in the cycle.Therefore, we used the spatial and temporal adaptive reflectance fusion model (STARFM) algorithm to fuse the data, supplemented the missing data, and obtained complete 2-day HJ-1 data for inversion.

STARFM Model
STARFM is a model based on reflectivity data [22].The main purpose is to use high spatial resolution, low temporal resolution reflectivity data, and low spatial resolution, high temporal resolution reflectivity data to obtain high spatial resolution, high temporal resolution reflectivity data.This algorithm was successfully applied in the seasonal variation map of vegetation, and good results were obtained [23][24][25].The algorithm has the following constructs:

STARFM Model
STARFM is a model based on reflectivity data [22].The main purpose is to use high spatial resolution, low temporal resolution reflectivity data, and low spatial resolution, high temporal resolution reflectivity data to obtain high spatial resolution, high temporal resolution reflectivity data.This algorithm was successfully applied in the seasonal variation map of vegetation, and good results were obtained [23][24][25].The algorithm has the following constructs: where L(x w/2 , y w/2 , T 0 ) is the HJ-1 predicted reflectivity value of time position; W is the size of the moving window, we used only pixels with spectral similarity to the central pixel for prediction in the window; M x i , y j , T k and L x i , y j , T k represent MODIS and HJ-1 reflectivity values of T k time in x i , y j window position, respectively; n represents the number of input pairs of images.
In the STARFM algorithm, the HJ-1 pixel size has been set as the basic unit, with the moving window size at 50, which is the least common multiple of the HJ-1 space resolution (30 m) and MODIS space resolution (500 m), in order to ensure that the pixels in the window are complete.The distance weight constant has been set at the half of the window size, which is a relatively balanced value.The uncertainty of HJ-1 and MODIS data has been set at 0.2% and 0.5%, which reflects the unreliable pixel ratio and is obtained from error statistics of pixel classification.The MODIS and HJ-1 data are used as the input images, and the weight values of neighboring relative center pixels are obtained by calculating the spectral, time and spatial distance weights; thereafter, the fusion image pixel values can be calculated.
A weighted average strategy is used in the STARFM algorithm.When the difference between the two sets of input images is too large, it may cause the "time smoothing" problem.The smaller the time interval, the lesser the difference between the input MODIS images, the more similar the optional spectral pixels, and the more accurate the forecast [26].Therefore, in the prediction, we selected the image pairs that were as close to the predicted image time as the prediction base images.

RTLSR Model
BRDF describes the direction of the surface reflection characteristics.In the bidirectional reflectivity models, a practical method called the kernel-driven model, which has a certain physical meaning of the linear combination of the core, is used to match the surface bidirectional reflection.RTLSR is the main algorithm for MODIS diversion and albedo products are currently providing albedo and bi-directional reflectivity of global scale to users.
The RTLSR algorithm is the weighted sum of two kernels, which are called "RossThick" and "LiSparse-Reciprocal".The RTLSR algorithm has the following constructs [27]: where K geo and K vol are geometrical optics kernel and bulk scattering kernel, which are functions of observing zenith angle θ i , solar zenith angle θ r and relative azimuth angle ϕ, respectively.f iso , f geo , and f vol are constant coefficients, which represent the proportion of isotropic scattering, geometric optical scattering, and bulk scattering.R(θ i , θ r , ϕ, λ) is the BRDF of the band.The RTLSR algorithm uses the linear regression to invert the best coefficients of the matching observation data, f iso , f geo , and f vol , then obtains the bi-directional reflectivity of random incident angle and observation angle by the extrapolation or interpolation of the kernel [28].
The bulk scattering kernel "RossThick" is suitable for describing dense homogeneous vegetation and leaf dip spherical distribution, proposed by Roujean [29]: where g is the phase angle defined by the following formula: The geometric optical kernel "LiSparseR" is suitable for describing a sparsely distributed canopy or other opaque geometry.It is a method based on the reciprocity principle, and it can partially improve the problems when extrapolating to the larger zenith angle without affecting the data matching ability.The expression is: where, where b is the vertical radius of the sphere; h is the horizontal radius of the sphere, and r is the height of the sphere.In the MODIS algorithm, h b = 2, b r = 1 [29].We used this algorithm, setting the HJ-1 fusion data to 16 days as a group and 8 days as a cycle, to match a set of kernel coefficients: f iso , f geo , and f vol .The least squares method was applied in matching.We set the solar zenith angle in the model to 45 • , which was consistent with the solar zenith angle of MOD43 NBAR product data, and the observation zenith angle to the zenith direction.According to the red and NIR kernel coefficients of each pixel in each cycle, and the relative azimuth ϕ, we obtained the azimuthal direction reflectivity HJ-1 NBAR, which was adjusted by the BRDF model for each pixel in the red and NIR bands.

SARIMA Model
Generally, historical data will have a strong relationship at the potential cycle time points, especially in economics.The SARIMA model is mainly used to identify the predictions of dependent variables, influenced by seasonal fluctuations and external events.The SARIMA model contains trend and seasonal variation, which has been widely used in different fields, such as economics, statistics, and RS, forecasting a certain parameter in time series with seasonality [30][31][32][33].
According to the definition, the SARIMA model is generally referred to SARIMA (p, d, q) × (P, D, Q) s .Where d and D are the order of the stepwise difference and the seasonal difference, respectively; p and q are the order of autoregressive and moving average, respectively; P and Q are the order of seasonal autoregressive and seasonal moving average, respectively; s is the seasonal period.The model can be expressed as follows: where, α t is the estimated residual at time t, a random process of Gaussian white noise.
The general matching process includes three steps: identification, estimation, and diagnostic checking.The crucial step in the SARIMA model is the choice of orders.The main methods used are autocorrelation functions (ACF) and partial autocorrelation functions (PACF), Akaike's information criterion (AIC), and so on.First, we used the ACF and PACF methods to select several possible orders.After modeling, the AIC criterion was used to select the best order model.The p value in the model was taken as 1, and the D value was taken as 1.
After the model was established, it needed to be tested.The method used here is the residual ACF test.In model testing, we are mainly concerned about whether the residual of the model is relevant, and normally distributed with a mean of zero.If the residuals of the SARIMA model are correlated and not normally distributed with a mean of zero, then the model can be further improved.On the contrary, the model fitting effect is good, and it can be considered that the model fully extracts the information of the sequence.
For the test of residual white noise, LB test was used to check whether the residual of the model is a white noise sequence.The test results are shown in Table 2, which indicates that when the residual sequence was delayed 1-12 orders, the P value of Q statistics is greater than 0.05.Therefore, at the significance level of 0.05, the original hypothesis is not rejected, meaning that the residual sequence is a white noise sequence, indicating that the fitted model had fully extracted the information in the time series.
For the residual distribution problem, we generated the model diagnosis report and analyzed the abnormal behavior of the research model, as shown in Figure 4. Figure 4a shows that the red KDE line follows the N (0, 1) line, which indicates that the residuals are normally distributed.Figure 4b shows that the residual error does not have autocorrelation, which indicates the model can help us understand the original time series data and predict the future values.
Using the HJ-1 NBAR data and MODIS NBAR data from interest areas 1 and 2 in 2011-2013 as the raw data, the model was used to predict the reflectivity of the red and NIR bands in 2014, and the results were compared with the surface measured data to discuss the prediction accuracy of the model.

Data Fusion Assessments
Before the data fusion process, we selected three images for the STARFM algorithm to illustrate the accuracy of the results and thus determine availability.The experimental data were HJ-1 data from areas of interest 1 and 2, and the corresponding time MOD09GA data.
Figure 5 shows three pairs of HJ-1 and MODIS images acquired on 26 April 2012, 4 May 2012, and 8 May 2012.The first row in Figure 5 shows three MODIS reflectivity images using band 2-1-4 as the red-green-blue composite, and the second row shows three HJ-1 reflectivity images using band 4-3-2 as the red-green-blue composite.The pairs of MODIS and HJ-1 images acquired on 26 April 2012 and 8 May 2012 were selected as base images, and the STARFM algorithm was used to predict the 4 May 2012 HJ-1 fusion image.
After STARFM algorithm calculation, we compared the predictions with the 4 May 2012 observed HJ-1 reflectivity image.Figure 6 shows the observed and predicted images produced by the STARFM algorithm in regions of interest 1 and 2. The boundary information of road, water, and vegetation can be expressed clearly in the fusion image, and there is high similarity between the predicted and the measured image.
There are many commonly used error evaluation indicators for model evaluation; we used several of these indicators to refine the results of the fusion evaluation, including root mean square error (RMSE), mean absolute percent error (MAPE), average absolute difference (AAD), and coefficient of determination (R 2 ).
RMSE was used to measure the deviation between the result and the reference value; the smaller the RMSE, the closer the result is to the reference value.
MAPE indicates the accuracy between the result and the reference value; if this value is small, it indicates high accuracy, according to the model predictive ability evaluation table proposed by Xia et al. [34].Generally, when MAPE ≤ 20, the forecast is valid.
Average absolute deviation (AAD) is the average of the absolute values of the difference between the result and the reference value.
R is a direct reflection of correlation between the result and the reference value.If R is positive, the correlation coefficient is positive, on the contrary, if R is negative, the two are negatively correlated.The larger the absolute value of R, the higher the degree of correlation, the greater the determination coefficient, the closer the relationship between the two sets of data.
Figure 7 shows scatter plots of the predicted and actual observations.Table 3 shows the analytical values for the predicted and actual observations.In the four scatter plots, the points are distributed on both sides of the y = x line.The distribution of the red and NIR band points of forest are closer to the straight line and more uniform, whereas the distributions of farmland are concentrated above the straight line.This means that the fusion values are slightly higher than the actual values.According to Table 3, the MAPE values of farmland are within 20%, and the R values are high, which indicates that the experimental results are effective.Compared to the reference image, the mean absolute difference is small, and the predicted value is highly correlated with the actual observed value.

Inversion Results
In the process of matching the kernel coefficient, the matching precision is higher when the data sampling direction is more uniform in the spatial distribution of the hemisphere.Figure 8a shows the red band solar zenith angle, observation zenith angle, and relative azimuthal distribution for a cycle (8 May to 22 May 2012) of HJ-1 data, selecting a 3 × 3 pixel window.The sampling direction of the data in the cycle is relatively uniform.To further determine the accuracy of the kernel coefficient matching, the kernel coefficient and the MOD43 BRDF kernel coefficient are used to create the direction reflectivity The results show that the fusion images have high spatial resolution, meanwhile, preserving periodic vegetation changes.The fusion images are close to the real observation images, and can be used to study the temporal variation of vegetation.Therefore, image fusion is performed to obtain complete 2-day HJ-1 reflectivity data.

Inversion Results
In the process of matching the kernel coefficient, the matching precision is higher when the data sampling direction is more uniform in the spatial distribution of the hemisphere.Figure 8a shows the red band solar zenith angle, observation zenith angle, and relative azimuthal distribution for a cycle (8 May to 22 May 2012) of HJ-1 data, selecting a 3 × 3 pixel window.The sampling direction of the data in the cycle is relatively uniform.
To further determine the accuracy of the kernel coefficient matching, the kernel coefficient and the MOD43 BRDF kernel coefficient are used to create the direction reflectivity distribution.The solar zenith angle is set to 45 • and the relative azimuth angle to the HJ-1 pixel relative azimuth angle.Figure 8b shows the red band inversion direction reflectivity for a cycle (8 May to 22 May 2012) of HJ-1 data and MOD43 BRDF data on 8 May 2012 in the farmland area of interest.The kernel coefficients of the matched nuclei and the MOD43 BRDF data of red band are 0.18559, 0, 0.0070382, and 0.209, 0.032, 0.049, respectively.The reflectivity of the inverted zenith direction is 0.1880 and 0.1832, respectively, whereas the MOD43 NBAR is 0.1828.Although the coefficients are different, the inversion of the reflectivity values is close.The inversion reflectivity of HJ-1 in the zenith direction is similar to that of MOD43 NBAR; it is proven that the inversion reflectivity result shows a small difference with the actual measurement result of MOD43.The results are within the error range.

Inversion Results
In the process of matching the kernel coefficient, the matching precision is higher when the data sampling direction is more uniform in the spatial distribution of the hemisphere.Figure 8a shows the red band solar zenith angle, observation zenith angle, and relative azimuthal distribution for a cycle (8 May to 22 May 2012) of HJ-1 data, selecting a 3 × 3 pixel window.The sampling direction of the data in the cycle is relatively uniform.To further determine the accuracy of the kernel coefficient matching, the kernel coefficient and the MOD43 BRDF kernel coefficient are used to create the direction reflectivity For the inversion results, i.e., the HJ-1 NBAR data, the MODIS NBAR product is used to compare with them, analyze their similarity and variance, and analyze the accuracy of the inversion results.
In the forest and farmland areas, the average reflectivity data of HJ-1 NBAR and MOD43 NBAR zenith direction were obtained.Figure 9 shows the reflectivity of HJ-1 NBAR and MODIS NBAR time curves.It can be observed that the average reflectivity of HJ-1 NBAR is similar to that of MODIS average reflectivity, and the fluctuation range of reflectivity is almost the same.However, it can also be seen that the average reflectivity of HJ-1 NBAR fluctuates more violently.In the forest area of interest, the reflectivity of the NIR band significantly increased during the growing period, but it does not maintain a high value smoothly.There is a small fluctuation during the period, and the reflectivity curve of the red band fluctuates vigorously, and there are obvious anomalies.In the farmland area of interest, the trend of reflectivity is similar to that of MODIS, but there are smaller fluctuations before the growth period.The wet and dry conditions of the bare soil may be the reason.This result indicates that the average reflectivity of HJ-1 NBAR is more suitable to reflect the subtle changes in the surface.
Similarly, we used the quantitative indicators to calculate the difference between the HJ-1 NBAR data and the measured values, and to analyze the accuracy of the inversion data.Figure 10 shows the scatter plots of the HJ-1 NBAR and MOD43 NBAR data.Table 4 shows the accuracy test result for the HJ-1 inversion data and MODIS inversion data.In the four scatter plots, the points are distributed on both sides of the y = x line.The MAPE values are within 20%, compared to the MOD43 NBAR, the mean absolute difference is small, and the predicted value is highly correlated with the actual observed value.
By comparison, it can be seen that the HJ-1 30 m resolution data are more sensitive than the MOD43 500 m resolution data.In the case of the same fluctuation trend, the HJ-1 reflectivity exhibits smaller fluctuations.Because of the underlying surface, there will be different surface reflectivities in the same vegetation growth conditions.Although there are more inflection points and a greater chance of outliers, the HJ-1 data can more realistically reflect the change in the surface during the test period.
As the results indicate, the combination of MODIS data and HJ-1 data is expected to be a method of improving the data quality of HJ-1 30 m NBAR time series.Certainly, compared with MODIS NBAR products, there is still a shortage of large-scale data usage in surface reflectivity inversion, when using the HJ-1 and MODIS joint data.However, with the use of small-scale data, the HJ-1 fusion data can highlight more detailed real surface conditions.

NBAR Time Series Analysis
Using the HJ-1 NBAR data, we can see the time series line graph, as shown in Figure 11.The reflectivity change in forestland is in accordance with the change in most of the vegetation reflectivity.At around the beginning of March each year, the vegetation starts to grow, and the NIR reflectivity significantly increases, remaining at around 0.25.Upon the arrival of winter, the vegetation withers and the NIR reflectivity decreases rapidly.The red band reflectivity is significantly lower than the NIR reflectivity, at around 0.05, and the fluctuation is small, increasing before spring and decreasing slowly during spring.As in farmland, the difference from most of the vegetation is due to human influence, as

NBAR Time Series Analysis
Using the HJ-1 NBAR data, we can see the time series line graph, as shown in Figure 11.The reflectivity change in forestland is in accordance with the change in most of the vegetation reflectivity.At around the beginning of March each year, the vegetation starts to grow, and the NIR reflectivity significantly increases, remaining at around 0.25.Upon the arrival of winter, the vegetation withers and the NIR reflectivity decreases rapidly.The red band reflectivity is significantly lower than the NIR reflectivity, at around 0.05, and the fluctuation is small, increasing before spring and decreasing slowly during spring.As in farmland, the difference from most of the vegetation is due to human influence, as

NBAR Time Series Analysis
Using the HJ-1 NBAR data, we can see the time series line graph, as shown in Figure 11.The reflectivity change in forestland is in accordance with the change in most of the vegetation reflectivity.At around the beginning of March each year, the vegetation starts to grow, and the NIR reflectivity significantly increases, remaining at around 0.25.Upon the arrival of winter, the vegetation withers and the NIR reflectivity decreases rapidly.The red band reflectivity is significantly lower than the NIR reflectivity, at around 0.05, and the fluctuation is small, increasing before spring and decreasing slowly during spring.As in farmland, the difference from most of the vegetation is due to human influence, as the surface is covered by bare soil before the growth period.The reflectivity changes in farmland are relatively stable.In the early stages of growth, the NIR band reflectivity increases, then decreases during the growth period, whereas the red band reflectivity decreases at the beginning of the growing season, and then increases during the growing season.In sum, the seasonal difference of forest and farmland albedo is mainly affected by the seasonal difference in plant growth.
This time series provided comparable dynamic data, avoiding the angular impact, to reflect the dynamic variation patterns of vegetated target areas.Therefore, we used the data to do the forecast analysis for disaster prediction, environmental control, etc.
We used the HJ-1 NBAR data and MODIS NBAR data of interest areas 1 and 2 as the raw data.For the HJ-1 NBAR data in the regions of interest, according to the average, we set a threshold of 5% and selected the pixels within the threshold to the new average, excluding mixed pixel effects.We implemented the SARIMA method using Eviews, which is freely available for data analysis, regression analysis, and forecasting.Observing the ACF and PACF figures, we selected the intervals of coefficients.After the seasonal difference and AIC screening, the model with the smallest AIC value was selected as the best, and the resulting model was determined in Table 5.
The white noise test of the model was effective.The residual sequences are purely random and the autocorrelation coefficients fall into the random interval, which indicates that the models are feasible and the prediction results are meaningful.The Shapiro normality test for the residuals produced p-values of 0.372 and 0.386, showing that the original hypothesis of following the normal distribution is rational, that is, the residuals exhibit a normal distribution.
For each land cover type and band, the predicted time series and the corresponding measured NBAR data are displayed in Figure 11.It can be seen from the figures that the HJ-1 predicted values and the measured values are very close.The calculated root mean square error is 8.36%, 7.58%, 8.87%, and 9.05%, indicating that the match is good.However, for the MODIS predicted values, especially in the forest, under the influence of mixed pixels, the predicted values are significantly higher than the measured values; there are large differences in small fluctuations.The HJ-1 NBAR predicted time series matches the measured NBAR values and reflects the seasonal changes in vegetation properties better; however, it is just an approximation of the MODIS NBAR predicted time series.In general, the predicted time series NBAR can represent the variation of interest area with a more reliable value compared with MODIS NBAR predicted data.The recursive estimation method is robust and reliable.
We used RMSE, MAPE, AAD, and R to evaluate the results, as shown in Table 6.The accuracy test values indicate that the HJ-1 predictions are closer to the real measured values than the MODIS predictions.Figure 12 shows scatter plots of the HJ-1 predicted and actual observations.In the scatter plots, the points are distributed on both sides of the y = x line.For the MODIS predicted values, the MAPE values are within 25% and the R value is 0.75.For the HJ-1 predicted values, the MAPE values are within 20%, the R value is 0.91, and the HJ-1 predicted value is highly correlated with the actual observed value.The HJ-1 predicted values match the measured data very well.Thus, the proposed HJ-1 inversion data and time series methods perform well for NBAR recursive estimation.
From these results, we found that (1) the fusion data had advantages of high spatial and temporal resolution, and matched HJ-1 data well; (2) the HJ-1 NBAR data adjusted by BRDF model were close to MOD43 data, and can reflect the annual variation characteristics of surface vegetation growth well; (3) the HJ-1 data had high-resolution features compared with the MOD43 NBAR data, reflecting the regional scale reflectivity spatial changes more accurately, and making related research and analysis more sophisticated; (4) compared to the forecast results of MODIS NBAR, the HJ-1 NBAR matched the surface measured data better, and can reflect the seasonal changes in vegetation properties better.

Conclusions
The inconsistent data quality in RS images, especially HJ-1 CCD data, is mainly influenced by environmental factors, and creates problems in the application of these data.In this paper, the mountain valley with farmland and forestland in North China is selected as the experimental area, and an improved method is presented to obtain the HJ-1 inversion NBAR data using linear matching of the RTLSR model, and then to predict reflectiv-

19 Figure 1 .
Figure 1.Location of the study area, Huailai County in northern Hebei Province of North China (a), and two interesting regions containing the measuring sites (b) in Huailai County.

Figure 1 .
Figure 1.Location of the study area, Huailai County in northern Hebei Province of North China (a), and two interesting regions containing the measuring sites (b) in Huailai County.

Figure 2 .
Figure 2. (a) Atmospheric correction contrast map, take the red band as an example; (b) spectral response function graph.Blue, green, red, and black curves in the figure correspond to blue, green, red, and near-infrared bands, respectively.

Figure 4 .
Figure 4. (a) Probability distribution diagram of model residual error, where KDE stands for kernel density estimation, N (0, 1) represents normal distribution, and Hist represents the model residual histogram; (b) autocorrelation function graph of residuals.

Figure 4 .
Figure 4. (a) Probability distribution diagram of model residual error, where KDE stands for kernel density estimation, N (0, 1) represents normal distribution, and Hist represents the model residual histogram; (b) autocorrelation function graph of residuals.

Figure 6 .
Figure 6.Observed image for 4 May ASD spectrometer 2012 and the predicted image: (a) observed image; (b) STARFM image.

Figure 6 .
Figure 6.Observed image for 4 May ASD spectrometer 2012 and the predicted image: (a) observed image; (b) STARFM image.

Figure 6 .
Figure 6.Observed image for 4 May ASD spectrometer 2012 and the predicted image: (a) observed image; (b) STARFM image.

Figure 7 .
Figure 7. Four scatter plots of fusion reflectivity and observed reflectivity (x, observed; y, fusion): (a) the distribution of the red band points of forest; (b) the distribution of the NIR band points of forest; (c) the distribution of the red band points of farmland; (d) the distribution of the NIR band points of farmland.

Figure 8 .
Figure 8.(a) Sun observation angle distribution, where the red circle, the green circle, and the pink circle are corresponding to the solar zenith angle, the original HJ-1 observation zenith angle, and the fusion HJ-1 data observation zenith angle, respectively; (b) directional reflectivity profile, where the red and blue lines are corresponding to the MOD43 and HJ-1 reflectivity, respectively.

Figure 7 .
Figure 7. Four scatter plots of fusion reflectivity and observed reflectivity (x, observed; y, fusion): (a) the distribution of the red band points of forest; (b) the distribution of the NIR band points of forest; (c) the distribution of the red band points of farmland; (d) the distribution of the NIR band points of farmland.

Figure 8 .
Figure 8.(a) Sun observation angle distribution, where the red circle, the green circle, and the pink circle are corresponding to the solar zenith angle, the original HJ-1 observation zenith angle, and the fusion HJ-1 data observation zenith angle, respectively; (b) directional reflectivity profile, where the red and blue lines are corresponding to the MOD43 and HJ-1 reflectivity, respectively.

Figure 8 .
Figure 8.(a) Sun observation angle distribution, where the red circle, the green circle, and the pink circle are corresponding to the solar zenith angle, the original HJ-1 observation zenith angle, and the fusion HJ-1 data observation zenith angle, respectively; (b) directional reflectivity profile, where the red and blue lines are corresponding to the MOD43 and HJ-1 reflectivity, respectively.

Figure 10 .
Figure 10.Four scatter plots of inversion reflectivity and observed reflectivity (x, observed MODIS NBAR; y, predicted HJ-1 NBAR): (a) the distribution of the red band points of forest; (b) the distribution of the NIR band points of forest; (c) the distribution of the red band points of farmland; (d) the distribution of the NIR band points of farmland.

Figure 10 .
Figure 10.Four scatter plots of inversion reflectivity and observed reflectivity (x, observed MODIS NBAR; y, predicted HJ-1 NBAR): (a) the distribution of the red band points of forest; (b) the distribution of the NIR band points of forest; (c) the distribution of the red band points of farmland; (d) the distribution of the NIR band points of farmland.

Figure 10 .
Figure 10.Four scatter plots of inversion reflectivity and observed reflectivity (x, observed MODIS NBAR; y, predicted HJ-1 NBAR): (a) the distribution of the red band points of forest; (b) the distribution of the NIR band points of forest; (c) the distribution of the red band points of farmland; (d) the distribution of the NIR band points of farmland.

Figure 11 .
Figure 11.Comparison between predicted and observed NBAR time series of the red and nearinfrared band data: (a) forest; (b) farmland.Each pair of the black, red, and blue curves in the subfigures correspond to red and near-infrared bands of the measured NBAR, HJ-1 predicted NBAR, and MODIS predicted NBAR, in which the red band curves have a relatively low value and the near-infrared curves have a relatively high value.

Figure 11 .
Figure 11.Comparison between predicted and observed NBAR time series of the red and nearinfrared band data: (a) forest; (b) farmland.Each pair of the black, red, and blue curves in the subfigures correspond to red and near-infrared bands of the measured NBAR, HJ-1 predicted NBAR, and MODIS predicted NBAR, in which the red band curves have a relatively low value and the near-infrared curves have a relatively high value.Appl.Sci.2022, 12, x FOR PEER REVIEW 17 of 19

Figure 12 .
Figure 12.Four scatter plots of HJ-1 predicted reflectivity and observed reflectivity (x, observed; y, predicted): (a) the distribution of the red band points of forest; (b) the distribution of the NIR band points of forest; (c) the distribution of the red band points of farmland; (d) the distribution of the NIR band points of farmland.

Figure 12 .
Figure 12.Four scatter plots of HJ-1 predicted reflectivity and observed reflectivity (x, observed; y, predicted): (a) the distribution of the red band points of forest; (b) the distribution of the NIR band points of forest; (c) the distribution of the red band points of farmland; (d) the distribution of the NIR band points of farmland.
40.19 • -40.21 • N and 115.78 • -115.80 • E) is covered by forest; there are four measuring sites in each region.
• N and 115.74 • -115.76 • E) is covered by farmlands, and corn is planted over 90% of the total area, and the south region (

Table 1 .
Experimental data and band comparison.

Table 2 .
White noise test results of fitting model residuals.

Table 3 .
Accuracy test of four data fusion results

Table 3 .
Accuracy test of four data fusion results

Table 3 .
Accuracy test of four data fusion results.

Table 4 .
Inversion result accuracy test.