Integration of Microwave and Optical/Infrared Derived Datasets from Multi-Satellite Products for Drought Monitoring

Drought is among the most common natural disasters in North China. In order to monitor the drought of the typically arid areas in North China, this study proposes an innovative multi-source remote sensing drought index called the improved Temperature–Vegetation–Soil Moisture Dryness Index (iTVMDI), which is based on passive microwave remote sensing data from the FengYun (FY)3B-Microwave Radiation Imager (MWRI) and optical and infrared data from the Moderate Resolution Imaging Spectroradiometer (MODIS), and takes the Shandong Province of China as the research area. The iTVMDI integrated the advantages of microwave and optical remote sensing data to improve the original Temperature–Vegetation–Soil Moisture Dryness Index (TVMDI) model, and was constructed based on the Modified Soil-Adjusted Vegetation Index (MSAVI), land surface temperature (LST), and downscaled soil moisture (SM) as the three-dimensional axes. The global land data assimilation system (GLDAS) SM, meteorological data and surface water were used to evaluate and verify the monitoring results. The results showed that iTVMDI had a higher negative correlation with GLDAS SM (R = −0.73) than TVMDI (R = −0.55). Additionally, the iTVMDI was well correlated with both precipitation and surface water, with mean correlation coefficients (R) of 0.65 and 0.81, respectively. Overall, the accuracy of drought estimation can be significantly improved by using multi-source satellite data to measure the required surface variables, and the iTVMDI is an effective method for monitoring the spatial and temporal variations of drought.


Introduction
Drought is usually a water deficit caused by an imbalance in the water supply and demand due to the lack of precipitation [1,2]. It can be commonly divided into three types: meteorological drought (water shortage caused by the imbalance between precipitation and evaporation amount for a long time), agricultural drought (insufficient soil moisture available to plants) and hydrological drought (the phenomenon when runoff does not reach the standard value or the water level of aquifer drops) [3,4]. Drought can cause serious social and economic impacts, so it is necessary to assess the variations in droughts [5]. An in-depth analysis of the spatial and temporal evolution characteristics of drought can provide a reference for integrated water resource management and drought response [6].
With the development of remote sensing technology, the amount of data is abundant, which can enable the procurement of a large range of land surface information [7,8]. Therefore, remote to determine the drought level to provide a real-time and reliable basis for future disaster prevention and mitigation work.

Study Area
Shandong Province is located in the eastern part of China and the downstream area of the Yellow River, with an area of 158,000 km 2 . It lies between E 114 • 09 -E 122 • 43 and N 34 • 22 -N 38 • 23 . The annual average temperature is between 11-14 • C, and the annual accumulated precipitation is 550-950 mm/year [43]. The research area is very sensitive to climate change due to the fragile ecosystem. In 2016, there was little rainfall in spring, and the precipitation recorded by 70 meteorological stations was the lowest in the same period of history. Crops suffered from water shortages caused by drought during the growing period, which has a seriously negative impact on economic development in the study area. The land types in Shandong Province are mainly divided into forestlands, grasslands, croplands, wetlands, built-up lands and other lands ( Figure 1). The main land cover types are built-up lands and croplands. According to the Shandong Statistical Yearbook, by the end of 2015, the land area for construction was 2.8201 million hectares and that of cultivated land was 7.611 million hectares. The total amount of water resources is 22.032 billion cubic meters, of which the total surface water resource is 12.118 billion cubic meters. The main lakes include Weishan Lake and Zhaoyang Lake. The rivers include the Yihe River and Majiahe River.
Water 2020, 12, x FOR PEER REVIEW 3 of 20 in order to determine the drought level to provide a real-time and reliable basis for future disaster prevention and mitigation work.

Study Area
Shandong Province is located in the eastern part of China and the downstream area of the Yellow River, with an area of 158,000 km 2 . It lies between E 114°09′-E 122°43′ and N 34°22′-N 38°23′. The annual average temperature is between 11-14 °C, and the annual accumulated precipitation is 550-950 mm/year [43]. The research area is very sensitive to climate change due to the fragile ecosystem. In 2016, there was little rainfall in spring, and the precipitation recorded by 70 meteorological stations was the lowest in the same period of history. Crops suffered from water shortages caused by drought during the growing period, which has a seriously negative impact on economic development in the study area. The land types in Shandong Province are mainly divided into forestlands, grasslands, croplands, wetlands, built-up lands and other lands ( Figure 1). The main land cover types are builtup lands and croplands. According to the Shandong Statistical Yearbook, by the end of 2015, the land area for construction was 2.8201 million hectares and that of cultivated land was 7.611 million hectares. The total amount of water resources is 22.032 billion cubic meters, of which the total surface water resource is 12.118 billion cubic meters. The main lakes include Weishan Lake and Zhaoyang Lake. The rivers include the Yihe River and Majiahe River.

Satellite Image Data and Data Preprocessing
The data used in this study include satellite remote sensing data (MODIS, FY-3B-MWRI), GLDAS SM, in situ meteorological data and agricultural statistical data. The dataset is specifically described, as shown below.

MODIS Data
The MODIS data (h27v05) from January to December 2016 used in this study include the MODIS LST product (MOD11A2, collection v006), surface reflectance product (MOD09Q1, collection v006), VI product (MOD13A2, collection v006) and albedo product (MCD43B3, collection v005) downloaded from the NASA Land Processes Distributed Active Archive Center (LP DAAC, https://lpdaac.usgs.gov/). The eight-day gridded level-3 product MOD09Q1 includes NIR and Red bands, with a spatial resolution of 250 m; the spatial resolution of MOD11A2, MOD13A2 and MCD43B3 is 1 km. The MODIS Reprojection Tool (MRT) was used for the projection transformation, resampling and format conversions of these products [44]. Because of the influence of sensor design, atmosphere and other factors when acquiring remote sensing images, there will be considerable noise in the image. Data will be partly missing as a result of clouds, so the acquired image data cannot accurately express the surface conditions. Therefore, it is necessary to reconstruct VI and LST using

Satellite Image Data and Data Preprocessing
The data used in this study include satellite remote sensing data (MODIS, FY-3B-MWRI), GLDAS SM, in situ meteorological data and agricultural statistical data. The dataset is specifically described, as shown below.

MODIS Data
The MODIS data (h27v05) from January to December 2016 used in this study include the MODIS LST product (MOD11A2, collection v006), surface reflectance product (MOD09Q1, collection v006), VI product (MOD13A2, collection v006) and albedo product (MCD43B3, collection v005) downloaded from the NASA Land Processes Distributed Active Archive Center (LP DAAC, https://lpdaac.usgs.gov/). The eight-day gridded level-3 product MOD09Q1 includes NIR and Red bands, with a spatial resolution of 250 m; the spatial resolution of MOD11A2, MOD13A2 and MCD43B3 is 1 km. The MODIS Reprojection Tool (MRT) was used for the projection transformation, resampling and format conversions of these products [44]. Because of the influence of sensor design, atmosphere and other factors when acquiring remote sensing images, there will be considerable noise in the image. Data will be partly missing as a result of clouds, so the acquired image data cannot accurately express the surface conditions. Therefore, it is necessary to reconstruct VI and LST using the Savitzky-Golay filtering method, which is a weighted average algorithm that uses smooth time series data as the moving window, and has been proven to more realistically reflect the growth and change in vegetation [45]. To eliminate deviations in the image space, a median filter that uses the median of neighboring pixels to replace outliers was implemented in the image space. Finally, the VI and LST were obtained for the monthly time scale.

FY3B Soil Moisture Product
Chinese FengYun (FY) series meteorological satellites were designed and manufactured to obtain meteorological observations, weather forecasts, and nautical aviation [46]. The MWRI aboard the FY-3B satellite is an instrument that obtains products by receiving electromagnetic radiation information inversion from both the horizontal and vertical directions for the Earth's atmosphere and the Earth's surface. The MWRI daily product of the global SM includes values retrieved by measurements from ascending (daytime) and descending (nighttime) half-orbits [47]. Through verification, it was found that the effective value proportion of ascending data was higher, so only the ascending data were used in this study. The data were downloaded from the National Satellite Meteorological Centre (http://nsmc.org.cn). The data projection is an Equal-Area Scalable Earth Grid (EASE-Grid) with a data format of HDF5 and a spatial resolution of 25 km [48].

Meteorological Data and Agricultural Statistical Data
The meteorological data, including temperature and precipitation data, were provided by the National Meteorological Information Center (http://data.cma.cn/). The monthly average temperature and monthly accumulated precipitation of nine national meteorological stations (Weifang, Zibo, Yanzhou, Huimin, Jinan, Dingtao, Laiwu, Feixian and Fushan) were selected, and the time range was from January to December 2016. The spatial distribution of nine national meteorological stations in Shandong Province is shown in Figure 1. Surface water, effective irrigation area and farmland area are from the Shandong Provincial Bureau of Statistics (http://www.stats-sd.gov.cn/).

GLDAS SM
The global land data assimilation system (GLDAS) is a global offline high-resolution ground simulation system, which drives a variety of offline surface models and combines ground observation data with satellite observation data to generate the optimal field of surface state and flux (https: //giovanni.gsfc.nasa.gov/). GLDAS-2.1 is one of the two components of the GLDAS version 2 (GLDAS-2) dataset. The dataset contains many variables, such as albedo, soil temperature, evapotranspiration and SM [49]. The SM content is divided into four grades: 0-10 cm, 10-40 cm, 40-100 cm and 100-200 cm. The use of GLDAS SM data with a depth of more than 10 cm may add some uncertainty to the assessment because FY3B-MWRI SM data can only obtain the land surface SM at 0-5 cm depth [50]. In addition, in many previous studies, the top 10 cm SM data from the surface models was used to verify various satellite SM data [50][51][52][53]. Therefore, 0-10 cm SM data from January to December were selected [54]. Figure 2 shows the flow chart for the TVMDI improvement, comparison and validation in this study. First, NDVI, LST and albedo with a resolution of 1 km in 2016 and SM data at 25 km were obtained from MODIS and FY3B-MWRI products. Then, we built a linear downscaling model and downscaled the SM data. iTVMDI and TVMDI were calculated using FY3B-MWRI-downscaled SM and MODSI data. In Section 4.1, GLDAS SM is based on a large number of observed data and multiple model simulations, which can represent the measured SM to some extent. Therefore, GLDAS SM was used to verify the accuracy of the results of the two models. We then used meteorological data Water 2020, 12, 1504 5 of 17 from nine meteorological stations in 2016 to verify the results of the two models and analyze the response capacity of the two models to meteorological drought (see Section 4.2). Then, we used surface water to test the monitoring capability of the two models for hydrological drought (see Section 4.3). Finally, we used the iTVMDI model to analyze the drought condition of Shandong Province in 2016 (see Section 4.4).

Improved Temperature-Vegetation-Soil Moisture Dryness Index (iTVMDI)
Amani et al. [34] proposed the concept of the TVMDI, which is based on the relationship among the vegetation, temperature and SM (see Appendix A for further details and calculation formula). The iTVMDI proposed in this study improves the TVMDI model. The TVMDI uses optical and infrared data to calculate the SM, which is largely affected by the data quality. Microwave remote sensing is considered to be the most effective way to obtain SM because of its strong physical mechanism [55]. Compared with the SM data calculated based on the NIR-Red spectral space, the microwave data are not only less affected by the weather, but also have better penetrability, the surface SM data under the vegetation can be obtained, and the higher temporal resolution makes the data more accurate. The FY3B-MWRI SM product uses a dual-channel inversion algorithm that can eliminate the influences of vegetation and roughness, and the inversion model of SM is based on the development of the parameterized microwave radiation model (multifrequency polarization (Qp) model). Therefore, SM derived from passive microwave remote sensing data was used to replace that derived from the optical and infrared data in the original model [56,57].
Since the spatial resolution of FY3B SM data is 25 km, it cannot be directly used to build a model with other data with a resolution of 1 km, so we downscaled the SM data to achieve the same scale for all data [58][59][60]. The downscaling method of multivariate statistical regression is currently one of the most commonly used downscaling methods with a simple calculation form and strong applicability, which can downscale passive microwave low-resolution SM data with the help of remote sensing data [61,62]. The principle of this method is to build a linear relationship between SM and its related variables at low spatial resolution and apply the relationship at a high spatial resolution to obtain high spatial resolution SM [63]. To reduce the dependence of LST, VI and albedo on the external environment and to facilitate comparisons among different regions, a linear model of SM and the normalized difference vegetation index (NDVI), LST and albedo was established, and the formula is as follows:

Improved Temperature-Vegetation-Soil Moisture Dryness Index (iTVMDI)
Amani et al. [34] proposed the concept of the TVMDI, which is based on the relationship among the vegetation, temperature and SM (see Appendix A for further details and calculation formula). The iTVMDI proposed in this study improves the TVMDI model. The TVMDI uses optical and infrared data to calculate the SM, which is largely affected by the data quality. Microwave remote sensing is considered to be the most effective way to obtain SM because of its strong physical mechanism [55]. Compared with the SM data calculated based on the NIR-Red spectral space, the microwave data are not only less affected by the weather, but also have better penetrability, the surface SM data under the vegetation can be obtained, and the higher temporal resolution makes the data more accurate. The FY3B-MWRI SM product uses a dual-channel inversion algorithm that can eliminate the influences of vegetation and roughness, and the inversion model of SM is based on the development of the parameterized microwave radiation model (multifrequency polarization (Qp) model). Therefore, SM derived from passive microwave remote sensing data was used to replace that derived from the optical and infrared data in the original model [56,57].
Since the spatial resolution of FY3B SM data is 25 km, it cannot be directly used to build a model with other data with a resolution of 1 km, so we downscaled the SM data to achieve the same scale for all data [58][59][60]. The downscaling method of multivariate statistical regression is currently one of the most commonly used downscaling methods with a simple calculation form and strong applicability, which can downscale passive microwave low-resolution SM data with the help of remote sensing data [61,62]. The principle of this method is to build a linear relationship between SM and its related variables at low spatial resolution and apply the relationship at a high spatial resolution to obtain high spatial resolution SM [63]. To reduce the dependence of LST, VI and albedo on the external environment and to facilitate comparisons among different regions, a linear model of SM and the normalized difference vegetation index (NDVI), LST and albedo was established, and the formula is as follows: S m = a 000 + a 001 A * + a 010 T * + a 100 NDVI * + a 002 A * 2 + a 020 T * 2 + a 200 NDVI * 2 + a 011 T * A * + a 101 A * NDVI * + a 110 NDVI * T * (1) where A * is the normalized surface albedo, NDVI * is the normalized VI, and T * is the normalized LST.
The process of normalization is as follows: Another improvement is to replace the PVI in the original model with the Modified Soil-Adjusted Vegetation Index (MSAVI). Because PVI is easily affected by the soil background, it is not the best VI for drought monitoring used in the TVMDI [64]. Some VIs, such as the Soil Adjusted Vegetation Index (SAVI) and MSAVI, have been developed to improve the sensitivity of vegetation by considering the influence of the atmosphere and soil effects [65]. By using MSAVI instead of PVI, the impact of the soil background can be effectively eliminated and the shortcomings of some vegetation indices that cannot correctly reflect the real situation of the land surface when the SM is high can be corrected. The MSAVI is calculated as follows: where R Red and R NIR are the reflectance values of the red and NIR bands, respectively. Therefore, the iTVMDI is constructed in a three-dimensional space integrating normalized LST, normalized downscaled SM and normalized MSAVI, and can be calculated using Equation (6).

Comparisons of Drought Indices and GLDAS SM
To evaluate the iTVMDI, the correlation coefficient (R) and p-value (P) between the iTVMDI and the GLDAS SM were selected. As shown in Table 1, the correlation coefficients between the SM and iTVMDI can exceed 0.62, with an average correlation coefficient of 0.73 and a significant correlation (P < 0.05). The slope of the fitting line is negative, indicating that they are negatively correlated. The average correlation coefficient between TVMDI and SM is 0.55, which cannot reach the level of the iTVMDI. Among all sites, Huimin, Fushan and Dingtao fail to pass the 0.05 significance test, and have the lowest correlation, with correlation coefficients of 0.49, 0.49 and 0.33, respectively. Therefore, compared with iTVMDI and TVMDI, iTVMDI is more stable than TVMDI. The reason for this is that the iTVMDI is based on microwave data, and the downscaled FY3B SM data are used to replace the NIR-Red spectral space data in the TVMDI model. Although the calculation of TVMDI is very simple and integrates the remote sensing data's optical and infrared bands, it is easily affected by land cover change, atmosphere and clouds. Therefore, the iTVMDI that is calculated by using the downscaling FY3B SM data is more reasonable. Furthermore, we calculated the monthly variation in GLDAS SM, iTVMDI and TVMDI, as shown in Figure 3. Because the drought index depends not only on SM, but also on VI and LST, the change in the drought index is not exactly the same as that of SM, but there should be a significant correlation between them. It is obvious that iTVMDI is more consistent with this feature than TVMDI.

Comparisons of Drought Indices and Meteorological Data
Drought often occurs in months with low precipitation and high temperatures [66]. In this study, the monthly precipitation data from nine weather stations were used to verify the ability of the iTVMDI to monitor meteorological drought. Figure 4 displays scatterplots of drought indices (iTVMDI and TVMDI) and rainfall. The correlation between the iTVMDI and precipitation is higher than that of the TVMDI, and the correlation coefficients R are 0.65 and 0.48 (p < 0.001), respectively. However, the iTVMDI and TVMDI are not strongly correlated with precipitation. The reason may be that the data selected in this study are on a monthly scale. After averaging different situations over a month, some information would be lost.

Comparisons of Drought Indices and Meteorological Data
Drought often occurs in months with low precipitation and high temperatures [66]. In this study, the monthly precipitation data from nine weather stations were used to verify the ability of the iTVMDI to monitor meteorological drought. Figure 4 displays scatterplots of drought indices (iTVMDI and TVMDI) and rainfall. The correlation between the iTVMDI and precipitation is higher than that of the TVMDI, and the correlation coefficients R are 0.65 and 0.48 (P < 0.001), respectively. However, the iTVMDI and TVMDI are not strongly correlated with precipitation. The reason may be that the data selected in this study are on a monthly scale. After averaging different situations over a month, some information would be lost.
Spring in Shandong usually lasts from March to May. Due to limited rainfall, drought often occurs in spring. Summer is from June to August, which is the monsoon season, and most of the rainfall in Shandong is concentrated in the summer monsoon period. Autumn starts in September and ends in November, and rainfall during this period is similar to that in spring. December to February of the next year is considered to be winter, with precipitation and temperature at their lowest levels in the whole year. The seasonal difference in precipitation in Shandong Province is large. Taking Jinan in Shandong Province as an example, the monthly cumulative precipitation, average temperature and iTVMDI were calculated ( Figure 5). The drought index showed a downward trend from winter to summer, and began to rise afterward. In summer, due to the high vegetation coverage, sufficient rainfall and the increase in SM, iTVMDI reached the lowest value at the end of summer (Figure 5b). After summer, the harvesting or withering of vegetation and the decrease in precipitation and SM led to the rapid recovery of the drought index to a high value state.   Water 2020, 12, x FOR PEER REVIEW 9 of 20 Spring in Shandong usually lasts from March to May. Due to limited rainfall, drought often occurs in spring. Summer is from June to August, which is the monsoon season, and most of the rainfall in Shandong is concentrated in the summer monsoon period. Autumn starts in September and ends in November, and rainfall during this period is similar to that in spring. December to February of the next year is considered to be winter, with precipitation and temperature at their lowest levels in the whole year. The seasonal difference in precipitation in Shandong Province is large. Taking Jinan in Shandong Province as an example, the monthly cumulative precipitation, average temperature and iTVMDI were calculated ( Figure 5). The drought index showed a downward trend from winter to summer, and began to rise afterward. In summer, due to the high vegetation coverage, sufficient rainfall and the increase in SM, iTVMDI reached the lowest value at the end of summer (Figure 5b). After summer, the harvesting or withering of vegetation and the decrease in precipitation and SM led to the rapid recovery of the drought index to a high value state.

Comparisons of Drought Indices with Surface Water Supply
The limited precipitation and SM lead to insufficient runoff, which causes hydrological drought. The increase in the surface water supply will increase the surface water. If the water supply is reduced, the regional temperature will increase, the air will dry, and the evaporation of water will accelerate as long as there is enough water in the soil, the ecological water demand will increase, and the water resources will be reduced [67]. At the same time, the iTVMDI will be increased [68,69]. Therefore, to evaluate the iTVMDI for hydrological drought monitoring, we selected the surface water supply per unit area of nine cities in Shandong Province from the Shandong Statistical Yearbook as auxiliary data and analyzed their correlation with the iTVMDI and TVMDI (Figure 6a). The results showed that the iTVMDI had a significant negative correlation with the surface water (R = −0.81, p < 0.01). With the increase in surface water, the iTVMDI will decrease. By comparing the correlation between TVMDI and surface water (R = −0.63), it can be found that iTVMDI is more sensitive to changes in surface water, which also means it is more accurate. Therefore, it can be concluded that iTVMDI is better than TVMDI in hydrological drought monitoring.

Comparisons of Drought Indices with Surface Water Supply
The limited precipitation and SM lead to insufficient runoff, which causes hydrological drought. The increase in the surface water supply will increase the surface water. If the water supply is reduced, the regional temperature will increase, the air will dry, and the evaporation of water will accelerate as long as there is enough water in the soil, the ecological water demand will increase, and the water resources will be reduced [67]. At the same time, the iTVMDI will be increased [68,69]. Therefore, to evaluate the iTVMDI for hydrological drought monitoring, we selected the surface water supply per unit area of nine cities in Shandong Province from the Shandong Statistical Yearbook as auxiliary data and analyzed their correlation with the iTVMDI and TVMDI (Figure 6a). The results showed that the iTVMDI had a significant negative correlation with the surface water (R = −0.81, P < 0.01). With the increase in surface water, the iTVMDI will decrease. By comparing the correlation between TVMDI and surface water (R = −0.63), it can be found that iTVMDI is more sensitive to changes in surface water, which also means it is more accurate. Therefore, it can be concluded that iTVMDI is better than TVMDI in hydrological drought monitoring.

Drought Monitoring Results in Shandong Province
It has been shown, over one year, that the iTVMDI has a relatively accurate drought monitoring capability, so the iTVMDI was used in this paper for drought monitoring in Shandong Province (Figure 7). The iTVMDI range of Shandong Province in 2016 was 0.24-0.85. The higher the value is, the drier the area is. By calculating the average iTVMDI value of each month in Shandong Province, it can be found that the iTVMDI was generally high in January, February and March, and the average iTVMDI of Shandong in those three months could reach 0.70. In particular, in March, the iTVMDI in some areas can reach the highest value of 0.85 in the whole year. According to the meteorological data, from February 14th to April 15th, the average precipitation across the whole province was only 5.0 mm, which was the lowest value in the same period since 1951. In addition, the water consumption of plant growth was large, which led to a rapid decrease in SM. In March, April and May, the high iTVMDI values were mainly distributed in the central and northern parts of Shandong Province, while Dezhou, Liaocheng, Heze, central and eastern Jining, south-central Linyi, and southern Zaozhuang were humid. In June, due to the harvest of winter wheat, the vegetation coverage decreased rapidly. Simultaneously, the temperature rose rapidly. Therefore, the iTVMDI in the northwestern, western and southwestern parts of Shandong Province, which are all plain areas, changed from humid to relatively dry. In July, August and September, it was relatively humid compared with the previous months, likely because of the abundant precipitation and the growth period of vegetation, and the average iTVMDI values were 0.53, 0.47 and 0.55, respectively.

Drought Monitoring Results in Shandong Province
It has been shown, over one year, that the iTVMDI has a relatively accurate drought monitoring capability, so the iTVMDI was used in this paper for drought monitoring in Shandong Province (Figure 7). The iTVMDI range of Shandong Province in 2016 was 0.24-0.85. The higher the value is, the drier the area is. By calculating the average iTVMDI value of each month in Shandong Province, it can be found that the iTVMDI was generally high in January, February and March, and the average iTVMDI of Shandong in those three months could reach 0.70. In particular, in March, the iTVMDI in some areas can reach the highest value of 0.85 in the whole year. According to the meteorological data, from February 14th to April 15th, the average precipitation across the whole province was only 5.0 mm, which was the lowest value in the same period since 1951. In addition, the water consumption of plant growth was large, which led to a rapid decrease in SM. In March, April and May, the high iTVMDI values were mainly distributed in the central and northern parts of Shandong Province, while Dezhou, Liaocheng, Heze, central and eastern Jining, south-central Linyi, and southern Zaozhuang were humid. In June, due to the harvest of winter wheat, the vegetation coverage decreased rapidly. Simultaneously, the temperature rose rapidly. Therefore, the iTVMDI in the northwestern, western and southwestern parts of Shandong Province, which are all plain areas, changed from humid to relatively dry. In July, August and September, it was relatively humid compared with the previous months, likely because of the abundant precipitation and the growth period of vegetation, and the average iTVMDI values were 0.53, 0.47 and 0.55, respectively. In order to further analyze the drought situation of different land types, we calculated the monthly change in the average iTVMDI of different land types over the whole year (Figure 8). For most of the year, the iTVMDI of croplands was the lowest among all land types. In addition, the iTVMDI of built-up lands and other lands was generally higher than that of other land types, which is also a feature that cannot be clearly shown in the previous drought index maps, and with the In order to further analyze the drought situation of different land types, we calculated the monthly change in the average iTVMDI of different land types over the whole year ( Figure 8). For most of the year, the iTVMDI of croplands was the lowest among all land types. In addition, the iTVMDI of built-up lands and other lands was generally higher than that of other land types, which is also a feature that cannot be clearly shown in the previous drought index maps, and with the increase in LST, the difference observed with other land types was more obvious. This phenomenon mainly occurs because, when the temperature rises, the decrease in surface water of built-up land and other land was accelerated, so compared with other land types, there was less surface water.

Discussion
Some studies have shown that an index that combines optical, infrared data with microwave data has a better drought monitoring capability [9,70]. Since the results of the TVMDI depend on the accuracy of the three indices, namely LST, VI and SM, we improved the TVMDI model using SM data derived from FY3B-MWRI to replace the MODIS SM. To further present differences between the iTVMDI and TVMDI, we used GLDAS SM to verify FY3B-MWRI SM and MODIS SM (Figure 9). It can be concluded that the correlation between FY3B-MWRI SM and GLDAS SM is higher than that between MODIS SM and GLDAS SM (p < 0.001), which also shows that the iTVMDI is better at drought monitoring in Shandong Province. The low correlation coefficient may be because the resolution of GLDAS SM (8 km) is inconsistent with the resolution of the two kinds of satellite data (1 km).

Discussion
Some studies have shown that an index that combines optical, infrared data with microwave data has a better drought monitoring capability [9,70]. Since the results of the TVMDI depend on the accuracy of the three indices, namely LST, VI and SM, we improved the TVMDI model using SM data derived from FY3B-MWRI to replace the MODIS SM. To further present differences between the iTVMDI and TVMDI, we used GLDAS SM to verify FY3B-MWRI SM and MODIS SM (Figure 9). It can be concluded that the correlation between FY3B-MWRI SM and GLDAS SM is higher than that between MODIS SM and GLDAS SM (P < 0.001), which also shows that the iTVMDI is better at drought monitoring in Shandong Province. The low correlation coefficient may be because the resolution of GLDAS SM (8 km) is inconsistent with the resolution of the two kinds of satellite data (1 km).  Figure 10 shows the SM in January, calculated using the MODIS NIR-Red band (Figure 10a), the downscaled FY3B SM (Figure 10b), the TVMDI (Figure 10c) and the iTVMDI (Figure 10d). To facilitate a comparison of the differences between the SM calculated using the NIR-Red band feature space, downscaled SM, TVMDI and iTVMDI in Shandong Province, all the calculation results were normalized. There are many differences between the SM data calculated by the MODIS NIR and Red bands and the FY3B-downscaled SM in terms of representing the drought situation in Shandong Province, China (Figure 10a, b). Compared with MODIS SM, the FY3B-downscaled SM can show the differences across different land covers. The TVMDI and iTVMDI show the same drought situations in Shandong Province overall. However, when the iTVMDI shows drought, the land surface actually shows normal or wet conditions because drought conditions are very sensitive to many factors, including the time lag effect between the monthly average iTVMDI and rainfall or surface variables, such as the LST and VI [70].  Figure 10 shows the SM in January, calculated using the MODIS NIR-Red band (Figure 10a), the downscaled FY3B SM (Figure 10b), the TVMDI (Figure 10c) and the iTVMDI (Figure 10d). To facilitate a comparison of the differences between the SM calculated using the NIR-Red band feature space, downscaled SM, TVMDI and iTVMDI in Shandong Province, all the calculation results were normalized. There are many differences between the SM data calculated by the MODIS NIR and Red bands and the FY3B-downscaled SM in terms of representing the drought situation in Shandong Province, China (Figure 10a,b). Compared with MODIS SM, the FY3B-downscaled SM can show the differences across different land covers. The TVMDI and iTVMDI show the same drought situations in Shandong Province overall. However, when the iTVMDI shows drought, the land surface actually shows normal or wet conditions because drought conditions are very sensitive to many factors, including the time lag effect between the monthly average iTVMDI and rainfall or surface variables, such as the LST and VI [70]. As a major province of grain production, most of Shandong is currently composed of croplands. We calculated the effective irrigated area of each city, which is the total area of paddy fields and irrigated land that are equipped with irrigation projects or equipment and can be irrigated normally, and this is an important index reflecting the drought resistance of cultivated land in China. According to the statistical yearbook data, the effective irrigated areas in Dezhou, Liaocheng, Heze and Jining were 492, 473, 631 and 473 thousand hectares, accounting for 60%, 68%, 66% and 61% of the local agricultural land area, respectively. The proportion of effective irrigated area in these four cities was the highest in the province. Because these areas were also the main production areas of winter wheat, the relatively complete irrigation facilities made the farmland water supply sufficient. The continuous water supply and intensive land use made these areas have relatively high average annual vegetation coverage (Figure 11a). However, it is still necessary to continue to build adequate irrigation measures in farmland areas, especially in Heze, Weifang and Qingdao, where the temperature is higher than in other areas (Figure 11b), in order to replenish groundwater and maintain SM in this area. Further research will continue to focus on drought risks to determine the impact of drought on agricultural productivity and livelihoods [20,71]. In combination with the variation relationship between drought index and measured surface water supply data, the relationship between the available amount of and demand for water resources is calculated in detail to improve the comprehensive utilization rate of water resources, which requires further exploration [72]. As a major province of grain production, most of Shandong is currently composed of croplands. We calculated the effective irrigated area of each city, which is the total area of paddy fields and irrigated land that are equipped with irrigation projects or equipment and can be irrigated normally, and this is an important index reflecting the drought resistance of cultivated land in China. According to the statistical yearbook data, the effective irrigated areas in Dezhou, Liaocheng, Heze and Jining were 492, 473, 631 and 473 thousand hectares, accounting for 60%, 68%, 66% and 61% of the local agricultural land area, respectively. The proportion of effective irrigated area in these four cities was the highest in the province. Because these areas were also the main production areas of winter wheat, the relatively complete irrigation facilities made the farmland water supply sufficient. The continuous water supply and intensive land use made these areas have relatively high average annual vegetation coverage (Figure 11a). However, it is still necessary to continue to build adequate irrigation measures in farmland areas, especially in Heze, Weifang and Qingdao, where the temperature is higher than in other areas (Figure 11b), in order to replenish groundwater and maintain SM in this area. Further research will continue to focus on drought risks to determine the impact of drought on agricultural productivity and livelihoods [20,71]. In combination with the variation relationship between drought index and measured surface water supply data, the relationship between the available amount of and demand for water resources is calculated in detail to improve the comprehensive utilization rate of water resources, which requires further exploration [72]. Water 2020, 12, x FOR PEER REVIEW 15 of 20

Conclusions
The TVMDI model is simple and efficient to calculate, but its results depend on SM, VI and LST obtained by optical and infrared remote sensing. Optical remote sensing is susceptible to weather and cannot obtain valid data in the presence of clouds. Microwave remote sensing is considered to be the most promising method of land surface SM remote sensing detection at present, given its short return period and high temporal resolution. In this paper, the original drought monitoring index TVMDI was improved by FY3B microwave remote sensing of SM. A novel drought index, iTVMDI, based on multi-satellite data fusion, was proposed in this study. Compared with the TVMDI model based on optical and infrared data, the iTVMDI model can improve the ability of drought monitoring by integrating optical/infrared data and microwave data. Furthermore, it also shows the potential of multi-source remote sensing data fusion in remote sensing and drought monitoring. In this study, iTVMDI was used to monitor the drought in Shandong Province in 2016, and GLDAS SM was used to verify the monitoring results. In addition, precipitation and surface water were used to assess the monitoring capacity of iTVMDI for meteorological and hydrological drought. The results proved that the iTVMDI model based on microwave data is more accurate, and the iTVMDI is better than TVMDI in meteorological drought and hydrological drought monitoring. Future research includes (1) improving the iTVMDI by using higher resolution microwave remote sensing SM data and (2) evaluating the ability of the iTVMDI for agricultural drought assessment.

Acknowledgments:
The MODIS data used in our study were originally downloaded from the NASA Land Processes Distributed Active Archive Center, located at the U.S. Geological Survey (USGS) Earth Resources Observation and Science (EROS) Center. The surface water data come from the Shandong Statistical Yearbook. The meteorological data come from the National Meteorological Information Center.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A
By increasing the SM, the vegetation coverage will increase. As a result, the VI value will rise, and the drought value will usually decrease. Similarly, when SM decreases and LST rises, vegetation will be affected by drought because there is not enough water in the soil. With the increase in LST

Conclusions
The TVMDI model is simple and efficient to calculate, but its results depend on SM, VI and LST obtained by optical and infrared remote sensing. Optical remote sensing is susceptible to weather and cannot obtain valid data in the presence of clouds. Microwave remote sensing is considered to be the most promising method of land surface SM remote sensing detection at present, given its short return period and high temporal resolution. In this paper, the original drought monitoring index TVMDI was improved by FY3B microwave remote sensing of SM. A novel drought index, iTVMDI, based on multi-satellite data fusion, was proposed in this study. Compared with the TVMDI model based on optical and infrared data, the iTVMDI model can improve the ability of drought monitoring by integrating optical/infrared data and microwave data. Furthermore, it also shows the potential of multi-source remote sensing data fusion in remote sensing and drought monitoring. In this study, iTVMDI was used to monitor the drought in Shandong Province in 2016, and GLDAS SM was used to verify the monitoring results. In addition, precipitation and surface water were used to assess the monitoring capacity of iTVMDI for meteorological and hydrological drought. The results proved that the iTVMDI model based on microwave data is more accurate, and the iTVMDI is better than TVMDI in meteorological drought and hydrological drought monitoring. Future research includes (1) improving the iTVMDI by using higher resolution microwave remote sensing SM data and (2) evaluating the ability of the iTVMDI for agricultural drought assessment. value of one. As shown in Figure A1, the line segment WD is the spatial representation line of TVMDI, and the point W indicates that the temperature in the region where the pixel is located is the lowest, the vegetation coverage is the highest, the SM content is the highest, and the drought value is the lowest. Notably, the value of the SM axis gradually decreases, because the lower the SM value is, the closer it is to the end with a TVMDI of one (Point D).
Water 2020, 12, x FOR PEER REVIEW 16 of 20 and the decrease in vegetation cover and SM, drought will be more serious and reach its maximum value of one. As shown in Figure A1, the line segment WD is the spatial representation line of TVMDI, and the point W indicates that the temperature in the region where the pixel is located is the lowest, the vegetation coverage is the highest, the SM content is the highest, and the drought value is the lowest. Notably, the value of the SM axis gradually decreases, because the lower the SM value is, the closer it is to the end with a TVMDI of one (Point D). Figure A1. The definition of the TVMDI. TVMDI is based on VI, LST, and SM as the axis of the threedimensional space. The TVMDI is defined as the square root of the square sum of VI, LST, and SM [34].
The TVMDI is defined as follows: where SM is calculated as follows: where Red R and NIR R are reflections in the red and NIR bands, respectively, after atmospheric correction, and M and b are the slope and intercept of the linear-fitting formula of the soil line in the NIR-Red feature space, respectively [24,73]. The PVI is defined as follows:  The TVMDI is defined as follows: where SM is calculated as follows: where R Red and R NIR are reflections in the red and NIR bands, respectively, after atmospheric correction, and M and b are the slope and intercept of the linear-fitting formula of the soil line in the NIR-Red feature space, respectively [24,73]. The PVI is defined as follows: where R Red and R NIR are reflections in the red and NIR bands, respectively. a and b refer to the slope and interception of the soil line equation, respectively. In this study, the values of each axis are normalized to the range from zero to √ 3 3 so that the drought index is in the range from zero to one, and the normalization formula is as follows: where P represents VI, SM or LST and P = (p 1 , p 2 , · · · , p n ). The maximum temperature is 349 K, and the minimum is 273 K because temperatures below 273 K are meaningless in drought research.