Development of a Multiple Linear Regression Model for Meteorological Drought Index Estimation Based on Landsat Satellite Imagery

Climate polarization due to global warming has increased the intensity of drought in some regions, and the need for drought estimation studies to help minimize damage is increasing. In this study, we constructed remote sensing and climate data for Boryeong, Chungcheongnam-do, Korea, and developed a model for drought index estimation by classifying data characteristics and applying multiple linear regression analysis. The drought indices estimated in this study include four types of standardized precipitation indices (SPI1, SPI3, SPI6, and SPI9) used as meteorological drought indices and calculated through cumulative precipitation. We then applied statistical analysis to the developed model and assessed its ability as a drought index estimation tool using remote sensing data. Our results showed that its adj.R2 value, achieved using cumulative precipitation for one month, was very low (approximately 0.003), while for the SPI3, SPI6, and SPI9 models, the adj.R2 values were significantly higher than the other models at 0.67, 0.64, and 0.56, respectively, when the same data were used.


Introduction
Climate change due to global warming has been noted worldwide and reported as causing new temperature and precipitation patterns [1], affecting agricultural water availability [2] and changing the effectiveness of irrigation systems [3]. Asia is also experiencing climate change due to global warming [4], and it has been predicted that precipitation patterns in China will include more frequent and severe precipitation patterns than those predicted in RCP8.5, a widely used climate change scenario model. Abnormal patterns have also been reported from Kazakhstan [5], where it was found that the increase in temperature due to global warming was a factor in the occurrence of flood damage. In addition, Weili et al. (2015) [6] analyzed changes in precipitation in Japan, from 1901 to 2012, and found that precipitation has decreased considerably in recent decades, while Mohammad et al. (2020) [7] reviewed annual climate change and evaluated associated trends in Iran from 1961 to 2010. Consequently, seasonal and regional variations in temperature and precipitation patterns appear to have increased, and the scale and frequency of damage due to droughts have also increased in some areas [8]. Drought occurs over a long period of time over a large area and is thus difficult to estimate. Drought is one of the most expensive natural disasters from which to temporally recover and is a potential risk to agriculture, water quality, and the economy. Additionally, the intensity of drought is

Study Area
This study used data from Boryeong, Chungcheongnam-do, South Korea, which has experienced frequent meteorological droughts over the past 10 years. Boryeong is a coastal city located in the mid-west of South Korea, and Boryeong Dam supplies water to Chungcheongnam-do. Since 2012, a countermeasure committee has been established at the central government level due to drought. In particular, the annual rainfall decreased to 1010 and 785 mm during 2014-2015 due to severe drought. The subject area and precipitation measuring points can be observed in Figure 1, while the trend line for SPI6, which was the standard index consulted for drought warnings in Korea during the study period, is shown in Figure 2.

Developing the MLR Model
MLR is a statistical technique that expresses the relationship between several independent variables and a dependent variable, representing the linear relationship as a single functional formula. The principle is the same as that of the simple linear regression, which reveals the relationship between one independent variable and a dependent variable; however, the dependent variable is generally affected by more than two independent variables in terms of explaining the most natural phenomena. MLR was applied because the accuracy of the regression model could be improved by selecting several independent variables [20]. The MLR model, using one dependent variable (y) and several independent variables (x ), uses the form shown in Equation (1) as follows:

Developing the MLR Model
MLR is a statistical technique that expresses the relationship between several independent variables and a dependent variable, representing the linear relationship as a single functional formula. The principle is the same as that of the simple linear regression, which reveals the relationship between one independent variable and a dependent variable; however, the dependent variable is generally affected by more than two independent variables in terms of explaining the most natural phenomena. MLR was applied because the accuracy of the regression model could be improved by selecting several independent variables [20]. The MLR model, using one dependent variable (y) and several independent variables (x i ), uses the form shown in Equation (1) as follows:

Developing the MLR Model
MLR is a statistical technique that expresses the relationship between several independent variables and a dependent variable, representing the linear relationship as a single functional formula. The principle is the same as that of the simple linear regression, which reveals the relationship between one independent variable and a dependent variable; however, the dependent variable is generally affected by more than two independent variables in terms of explaining the most natural phenomena. MLR was applied because the accuracy of the regression model could be improved by selecting several independent variables [20]. The MLR model, using one dependent variable (y) and several independent variables (x i ), uses the form shown in Equation (1) as follows: where β i(i=1,2,···k) represents a regression coefficient for the independent variable x i(i=1,2,···k) , y refers to the dependent variable, and C is the constant of the regression equation. The ordinary least squares method, which is commonly used when estimating regression coefficients, was used, and all independent variables used in the analysis were included in the regression equation by applying the simultaneous input method.

Model Assessment and Selection
The models developed by applying MLR analysis were evaluated in two ways. The first method measured the error of the SPI derived using each model in comparison with the actual SPI, while the second method compared the coefficient of determination of the derived regression equation. For the error evaluation index, the root mean squared error (RMSE) and mean absolute error (MAE) were calculated and compared. It has the characteristic of returning the error in a unit similar to that of the actual value. Regarding the RMSE, the square mean of the residuals of the estimated and actual values was calculated and square-rooted. RMSE is the most commonly used error evaluation index (Equation (2)). Compared to MAE, it has the characteristic of being sensitive to models with larger error values. MAE is the mean of the difference between actual and estimated values converted to an absolute value. It has an advantage over the RMSE when analyzing data with several outliers (Equation (3)).
where m indicates the number of dataset days used for testing, and y i stands for the estimated value of the drought index for the applicable date obtained from the developed regression model.ŷ i denotes the actual drought index for that date. The coefficient of determination, which evaluates the performance of the derived regression model, and confirms its reliability, verifies the R 2 and calculates adj.R 2 (Equation (4)) values.
Here, SSE and SST are the sum of squares of error and total sum of squares of each model, n is the number of data, and k is the number of independent variable types R 2 is an index that measures the degree to which the estimated linear model is suitable for the given data. It is generally interpreted as the explanatory power of the model, but increases as independent variables are added [21]; thus, when using MLR, R 2 and adj.R 2 were confirmed together. The 24 developed regression models were classified based on the SPI used as the dependent variable, and the performance of the models was compared. Figure 3 illustrates an overall schematic diagram of the research method. the number of data, and k is the number of independent variable types R² is an index that measures the degree to which the estimated linear model is suitable for the given data. It is generally interpreted as the explanatory power of the model, but increases as independent variables are added [21]; thus, when using MLR, R² and adj.R² were confirmed together. The 24 developed regression models were classified based on the SPI used as the dependent variable, and the performance of the models was compared. Figure 3 illustrates an overall schematic diagram of the research method.

Data
Drought index, satellite image data, and meteorological data were collected from reliable institutions to build a dataset to be used for the development of a MLR model. Drought index and meteorological data were collected from the Korea Meteorological Administration; satellite image data, i.e., Landsat 5 and Landsat 8 satellite images, were collected from the USGS. Subsequently, through QGIS, mean values for the three indices (NDVI, NDMI, and land surface temperature (LST)) used in determining drought were calculated. The raw data collated for this study have been itemized in Table 1. The study period was from June 2010 to September 2019, when droughts occurred frequently in Boryeong, and data were classified for a total of 76 days, facilitating their securement and development. We developed 24 models for this study by applying MLR analysis to each dataset and classifying them in groups of three, as shown in Figure 4. The SPI was used as the first classification criterion and was assigned as a dependent variable. Since four SPI types were used in this study (SPI1, SPI3, SPI6, and SPI9), there was a total of four cases.

Drought Indices
A drought index can express drought damage severity quantitatively, and several types have been developed to date. We used the drought index SPI as a dependent variable for regression model development. This drought index was developed using the idea that drought starts from a lack of precipitation [22]. It is the drought index most widely used to indicate drought severity. SPEI is a drought index that uses precipitation and evapotranspiration [23]; the method is similar to that of SPI, but is calculated by excluding cumulative evaporation from cumulative precipitation. Palmer expressed the depth of drought as a function of water shortage and the water shortage period via the Palmer Drought Severity Index (PDSI) [24]. Onyutha developed standardized non-parametric indices of precipitation and evaporation (SNIPE) The weakness of the drought index was compensated [25]. The drought index used as the dependent variable for regression model development is the Standardized Precipitation Index (SPI). In this study, four SPI types (SPI1, SPI3, SPI6, and SPI9) were collected, with the index number denoting the number of months (30 days per month) of cumulative precipitation used in its development. South Korea prepares its four weather drought warning stages using SPI6: mild drought (SPI6 < −1.0), moderate drought (<−1.5), severe drought (<−2.0), and extreme drought (<−2.0, lasting for >20 days). SPI is an index that expresses the optimal fit of precipitation to a probability distribution, and is one of the most widely used drought indices in modern times as suggested by Mckee et al. SPI is calculated using a standardized value of cumulative precipitation over a given period based on a 30 y precipitation record [26]. Because it has a variety of time scales, The second classification criterion involved using the characteristics of the dataset, and this was applied as an independent variable. In this study, climate data, remote sensing data, or both were used as independent variables (all-type data). The study was conducted by dividing it into three data types, with the final classification being based on absolute SPI values, with cases divided into those with an absolute SPI value > 1 and those where it was < 1. We applied this distinction as the number of days when the SPI did not exceed 1 involved more than half of the built datasets. Moreover, in the evaluation stage, the performance of the regression model developed using undifferentiated data was rated as poor, owing to the data from dates with absolute SPI values > 1.
From data constructed in such a manner that the performance of the model could be evaluated more objectively, a model was developed by subjecting approximately 80% of the target period data to MLR analysis. For the rest, the performance was assessed by using values estimated through the model and comparing them with actual values. The number of days used for training each data classification and the testing dataset days are listed in Table 2.

Drought Indices
A drought index can express drought damage severity quantitatively, and several types have been developed to date. We used the drought index SPI as a dependent variable for regression model development. This drought index was developed using the idea that drought starts from a lack of precipitation [22]. It is the drought index most widely used to indicate drought severity. SPEI is a drought index that uses precipitation and evapotranspiration [23]; the method is similar to that of SPI, but is calculated by excluding cumulative evaporation from cumulative precipitation. Palmer expressed the depth of drought as a function of water shortage and the water shortage period via the Palmer Drought Severity Index (PDSI) [24]. Onyutha developed standardized non-parametric indices of precipitation and evaporation (SNIPE) The weakness of the drought index was compensated [25]. The drought index used as the dependent variable for regression model development is the Standardized Precipitation Index (SPI). In this study, four SPI types (SPI1, SPI3, SPI6, and SPI9) were collected, with the index number denoting the number of months (30 days per month) of cumulative precipitation used in its development. South Korea prepares its four weather drought warning stages using SPI6: mild drought (SPI6 < −1.0), moderate drought (<−1.5), severe drought (<−2.0), and extreme drought (<−2.0, lasting for >20 days). SPI is an index that expresses the optimal fit of precipitation to a probability distribution, and is one of the most widely used drought indices in modern times as suggested by Mckee et al. SPI is calculated using a standardized value of cumulative precipitation over a given period based on a 30 y precipitation record [26]. Because it has a variety of time scales, SPI is also used for drought monitoring, early warning, and drought severity estimation. Hydrological drought monitoring is possible with the use of a long-term SPI [27].

Climate Data
The data used as independent variables in this study were divided into two categories. The first type, point data measured at weather stations on the ground, were collected from two measuring platforms: automatic weather systems (AWSs) and automated synoptic observing systems (ASOSs). A total of 77 datasets was secured from June 2010 to June 2019 based on the date on which remote sensing data were collected. Average wind speed and daily precipitation data were collected from the AWS, which is an observation system involved in preventing natural disasters caused by weather phenomena such as typhoons, floods, and droughts. Based on the location of the AWS, a Thiessen polygon was created and replaced with the area data. ASOS is a ground observation system that is carried out simultaneously at all stations to determine the atmospheric conditions at a set time. It was used to collect certain weather elements that were not observed by the AWS, and the ASOS within Boryeong was selected. The ASOS data collected in this study were the local atmospheric pressure, average relative humidity, and average time of sunshine.

Remote Sensing Data
Remote sensing data refer to data collected remotely through satellites, and in this study, Landsat satellite imagery was used. The Landsat series of satellites supply photographic imagery covering the entire earth. These satellites were jointly developed by NASA and USGS, and eight satellites (Landsat 1 (1972) to Landsat 8 (2013)) have been launched in this series so far. Landsat satellite data are characterized by high quality and easy acquisition. They provide data in bands covering various wavelengths (see Table 3) and the required index can be calculated using this. For the collected satellite images, data on the days without cloud cover over Boryeong City were used from June 2010 to June 2019. The band types provided by each satellite can be verified through Table 3. The index used as an independent variable can be calculated using the Landsat bands required for NDVI, NDMI, and LST. Figure 5 shows the remote sensing data of Boryeong on 13 June 2019, which were obtained from Landsat satellite images. In this study, the three area-averaged indicators-NDVI, NDMI, and LST-were computed using Landsat 5 and Landsat 8 satellite imagery. phenomena such as typhoons, floods, and droughts. Based on the location of the AWS, a Thiessen polygon was created and replaced with the area data. ASOS is a ground observation system that is carried out simultaneously at all stations to determine the atmospheric conditions at a set time. It was used to collect certain weather elements that were not observed by the AWS, and the ASOS within Boryeong was selected. The ASOS data collected in this study were the local atmospheric pressure, average relative humidity, and average time of sunshine.

Remote Sensing Data
Remote sensing data refer to data collected remotely through satellites, and in this study, Landsat satellite imagery was used. The Landsat series of satellites supply photographic imagery covering the entire earth. These satellites were jointly developed by NASA and USGS, and eight satellites (Landsat 1 (1972) to Landsat 8 (2013)) have been launched in this series so far. Landsat satellite data are characterized by high quality and easy acquisition. They provide data in bands covering various wavelengths (see Table 3) and the required index can be calculated using this. For the collected satellite images, data on the days without cloud cover over Boryeong City were used from June 2010 to June 2019. The band types provided by each satellite can be verified through Table  3. The index used as an independent variable can be calculated using the Landsat bands required for NDVI, NDMI, and LST. Figure 5 shows the remote sensing data of Boryeong on 13 June 2019, which were obtained from Landsat satellite images. In this study, the three area-averaged indicators-NDVI, NDMI, and LST-were computed using Landsat 5 and Landsat 8 satellite imagery.   NDVI analyzes the difference between the reflectance at the near-infrared (NIR) and red wavelengths and is the most widely used vegetation-related index. In a healthy vegetation area, the red wavelength is absorbed and the near-infrared wavelength has a high reflectance. Conversely, in the case of soil without vegetation, the reflectance in the red region is high, but that in the near-infrared region is low [28]. To emphasize this characteristic, the NDVI-which ranges from 1 to −1, with a higher value indicating healthy vegetation-is denoted as shown in Equation (5) NDMI is used to determine vegetation moisture content. It focuses on removing changes due to the leaf's internal structure and dry matter content in the vegetated area and explores vegetation moisture content by highlighting the difference between NIR and short-wavelength infrared (SWIR) measurements. The reflectance of SWIR is inversely proportional to the moisture content of the leaf, and the NDMI is represented as shown in Equation (6) Landsat provides the amount of energy observed for each channel, which is used to calculate LST data for that digit number (DN). For LST calculation, the data were converted into the actual amount of radiation through an equation provided by the USGS [29]. Landsat 5's LST was calculated by Water 2020, 12, 3393 9 of 16 substituting Equation (7) for Band 6, and for Landsat 8, LST was calculated by substituting Equation (8) for Bands 10 and 11.
where L λ is the spectral radiation amount reaching the sensor, Q cal shows the DN of the pixel unit analyzed in the image data, L Minλ represents the spectral radiation amount when Q cal is zero, and L Maxλ denotes the spectral radiation amount when Q cal = Q calmax . Q calmax and Q Calmin are the values expressed in DN units after quantifying the minimum and maximum radiation amounts, respectively.
where L λ represents the amount of spectral radiation reaching the sensor, and M L is the radiance multiplicative scaling factor for the band. Q cal represents the DN value of the pixel, and A L denotes the radiance additive scaling factor for the band. The radiation calculation was used to determine the "brightness" temperature, as shown in Equation (9) below: where T indicates the brightness temperature (K), and coefficients K 1 and K 2 (as Watts/(m 2 ·srad·µm)) represent correction factors provided by the USGS, as shown in Table 4. In this study, Landsat 8 Band 10 data were used for LST calculations. We also needed to calculate the emission rate (ε), which was determined as shown in Table 5 using the NDVI range. The USGS recommends not relying on values calculated using Landsat 8 Band 11 for LST calculation due to its higher levels of uncertainty; therefore, in this study, LST was calculated using data from Band 10 only [30]. Finally, LST values were calculated using a plugin provided by QGIS, which applied Equation (10) provided below (unit: K): The construction process applied to create each dataset has been illustrated in the schematic shown in Figure 6.
Finally, LST values were calculated using a plugin provided by QGIS, which applied Equation (10) provided below (unit: K): The construction process applied to create each dataset has been illustrated in the schematic shown in Figure 6.

Coefficient of Determination
The SPI value was estimated using each of the constructed MLR models. Table 6 presents the coefficients of the regression model developed by applying the previously suggested method to the all-type dataset. The t-value indicates the significance of each coefficient; the larger the absolute value, the greater the significance. The overall model summary, where R² and adj.R² values were verified, is presented in Table 7 and Each case with the highest coefficient of determination of the drought index was shaded. The F-value indicates the significance of the regression equation, and the

Coefficient of Determination
The SPI value was estimated using each of the constructed MLR models. Table 6 presents the coefficients of the regression model developed by applying the previously suggested method to the all-type dataset. The t-value indicates the significance of each coefficient; the larger the absolute value, the greater the significance. The overall model summary, where R 2 and adj.R 2 values were verified, is presented in Table 7 and Each case with the highest coefficient of determination of the drought index was shaded. The F-value indicates the significance of the regression equation, and the larger the value, the greater the significance. The models with |ŷ| > 1 exhibited higher adj.R 2 values than those with |ŷ | < 1 in general, and datasets with more variables had more significance.
However, there were exceptions. First, SPI1 did not markedly vary across most models. The adj.R 2 value of the |ŷ| > 1 ( = −0.00062) model was low at −0.19763. Moreover, SPI6 of the |ŷ| > 1 model (the remote sensing dataset), which used the least independent variables, had the highest adj.R 2 value at 0.64204. In the |ŷ| > 1 model, the adj.R 2 value for SPI3, which used all the variables, was the highest at 0.695654. The |ŷ| > 1 for the same drought index and climate dataset model had the second-highest value at 0.672715; however, the adj.R 2 value for the |ŷ| < 1 model was low, and this poor performance indicated that it was not reasonable to use a regression model that achieved such results.  RMSEs and MAEs were calculated and compared to evaluate the drought index estimation ability of the models developed in this study. As the R 2 and adj.R 2 values of |SPI| < 1 were not satisfactory, a model with different criteria needed to be selected to present the regression equation in this case. The RMSE and MAE, which are used as error indicators, were applied to calculate the comparative residuals between the estimated and actual SPI values calculated using the developed model. Table 8 shows the RMSE and MAE of each model, and the error evaluation indicators with the best performance were shaded. All the error indicators had values lower than the |ŷ| < 1 dataset, which indicated that using |ŷ| < 1 would result in better estimations. The |ŷ| < 1 climate dataset demonstrated better estimation ability than other dataset types, and SPI6 had the lowest RMSE and MAE of all the indices. Furthermore, the |ŷ| < 1 dataset performed better in the SPI6 and SPI9 models, which have been used to determine drought using relatively long-term rainfall data, than it did in the SPI1 and SPI3 models, which use short-term precipitation data. In fact, SPI9 had the highest estimation ability of the indices using the |ŷ| < 1 climate dataset, followed by SPI1, SPI6, and SPI3. However, notably, the R 2 and adj.R 2 values of the SPI1 model were unsatisfactory.

Best Model Selection
The size of adj.R 2 was selected as a criterion, based on the results shown in Table 8, for the |ŷ| > 1 models, while for the |ŷ| < 1 models, RMSEs were compared, with low values being preferred. A graph comparing the adj.R 2 and RMSE values calculated for each model is shown in Figure 7, with the regression model for each SPI considered in this study presented in Table 9 and notable results were shaded.
Water 2020, 12, x FOR PEER REVIEW 12 of 15

Best Model Selection
The size of adj.R² was selected as a criterion, based on the results shown in Table 8, for the |y | > 1 models, while for the |y | < 1 models, RMSEs were compared, with low values being preferred. A graph comparing the adj.R² and RMSE values calculated for each model is shown in Figure 7, with the regression model for each SPI considered in this study presented in Table 9 and notable results were shaded.
These results showed that all |y | < 1 models had their lowest RMSE when only the climate dataset was used and that the remote sensing data did not significantly influence their SPI estimate results. With respect to the |y | > 1 models selected using adj.R² values, all models that used the remote sensing data performed satisfactorily, except for SPI1. Thus, models that used only the remote sensing data were selected for SPI6. Our results showed that the SPI3 and SPI6 models had high adj.R² values compared to other models. However, their RMSEs were lower than those of the other models, and the RMSE calculated for the SPI3 | y | > 1 model being the highest at 4.65. SPI9 had the best performance when estimation ability of the model was assessed using RMSEs. In addition, unlike other models with a value >1, the SPI9 |y | > 1 model had high adj.R² and low RMSE values. However, the fact that the adj.R² value of the SPI9 |y |>1 model was the lowest among the selected models had to be taken into account.    These results showed that all |ŷ| < 1 models had their lowest RMSE when only the climate dataset was used and that the remote sensing data did not significantly influence their SPI estimate results. With respect to the |ŷ| > 1 models selected using adj.R 2 values, all models that used the remote sensing data performed satisfactorily, except for SPI1. Thus, models that used only the remote sensing data were selected for SPI6. Our results showed that the SPI3 and SPI6 models had high adj.R 2 values compared to other models. However, their RMSEs were lower than those of the other models, and the RMSE calculated for the SPI3 |ŷ| > 1 model being the highest at 4.65. SPI9 had the best performance when estimation ability of the model was assessed using RMSEs. In addition, unlike other models with a value >1, the SPI9 |ŷ| > 1 model had high adj.R 2 and low RMSE values. However, the fact that the adj.R 2 value of the SPI9 |ŷ|>1 model was the lowest among the selected models had to be taken into account.

Discussion
This study developed an MLR model and evaluated the possibility of using Landsat remote sensing data to compensate for the weakness of the existing research that estimates only a very small range of drought. To evaluate the applicability of remote sensing data for drought estimation, we extracted data to be used in drought indices from Landsat satellite imagery covering Boryeong, and developed a model for the application of MLR analyses. Most of the |ŷ| > 1 model results, excluding SPI1, yielded a higher adj.R 2 than the models using remote sensing and climate data. Thus, remote sensing data are expected to be effective when estimating the lack of accumulated precipitation for 3-6 months. However, the results of this study revealed limitations. First, among the models using remote sensing data, the models with SPI1 as the dependent variable showed poor overall performance. NDVI and NDMI included in the remote sensing dataset type are affected by vegetation and soil moisture, respectively, and these indices are also affected by long-term meteorological phenomena [31].
We found that the differences between the coefficients of determination and error indices were large and related to the range of |ŷ|. For the constructed datasets, Landsat satellite imagery data were collected for cloudless days over Boryeong, rather than at regular intervals. This approach-in addition to the fact that a Landsat satellite only observes the same area for approximately 16 days [32]-made building data using regular intervals difficult. [33] Consequently, accounting for the weather between dates established using only the SPI and remote sensing data was difficult. The use of TERRA/AQUA MODIS data-which are updated on a daily basis-or GEO-KOMPSAT satellite imagery-which continually observes only Korea-would have achieved satisfactory performance regardless of the |ŷ| range. Based on this research methodology, a more diverse model should be developed and applied to other regions of similar size. It will then be possible to expand the research results by applying a more versatile model to national (and continental) scale units. Notably, SPI is calculated using cumulative precipitation, although NDVI and NDMI are more related to soil moisture, suggesting the need for a follow-up study that examines a drought index other than the SPI as an independent variable. This follow-up study may apply indices selected based on their ability to estimate agricultural or hydrological drought, moving beyond simple meteorological drought estimation. Lastly, the results of this study suggest that the application of deep learning (rather than machine learning) techniques, such as MLR analysis, can improve performance and nonlinear expressiveness, resulting in the development of more accurate and reliable models. [34]. Based on this research methodology, a more diverse model should be developed and applied to other regions of a similar size to expand the research results by developing a highly versatile model through application to entire countries and continents.

Conclusions
Water shortages can be addressed and drought-related damage minimized by estimating drought and establishing countermeasures in advance. An important task in water resource management is preparing for extreme drought in the near future by exploring satellite remote sensing data that are suitable for drought estimation research, and creating methodology that can develop better models based on research results.