Managing Agricultural Water Considering Water Allocation Priority Based on Remote Sensing Data

Abstract

To fairly distribute limited irrigation water resources in arid regions, a water allocation priority evaluation method based on remote sensing data was proposed and integrated with an optimization model. First, the water supply response unit was divided according to canal system conditions. Then, a spatialization method was used for generating spatial agricultural output value (income from planting industry) and grain yield (yield of food crops) with the help of NDVI and the potential yield of farmland. Third, the AHP-TOPSIS method was employed to calculate the water allocation priority based on the above information. Finally, the evaluation results were integrated with a nonlinear multiobjective model to optimally allocate agricultural land and water resources, considering the combined objective of minimum envy and proportional fairness. The method was applied to Hetao irrigation area, an arid agriculture-dominant region in Northwest China. After solving the model, optimization alternatives were obtained, which indicate that: (1) the spatial method of agricultural output value can improve the accuracy by around 16% compared with the traditional method, and the spatial method of grain yield also have good accuracy (MAPE = 14.66%); (2) the rank of water allocation priority can reflect more spatial information, and provide practical decision support for the distribution of water resources; (3) the envy index can better improve the efficiency of an allocation system compared to the Gini coefficient method.

1. Introduction

Fresh water is an essential resource. It not only helps to maintain ecological health, but also promotes social development [1]. Agriculture is the largest water use sector, accounting for approximately 70% of freshwater use globally per year [2]. The ever-increasing food demand and decreased availability of water agriculture, calls for effective agricultural water management [3]. Especially in arid irrigation districts, when water scarcity occurs, water managers are required to make difficult decisions about what land has higher priority and how to effectively allocate limited water resources [4].

Priority ranking among different water users has been proven to be a useful way to help water managers identify the importance of water users [4,5]. A series of priority evaluation models for different water use sectors such as agriculture and industry, considering various factors such as social economy, ecological environment, and policies, have been developed [6,7,8]. For example, RazaviToosi and Samani presented an assessment model including social, economic, managerial and environmental clusters to evaluate the priority of five watersheds in Iran [9]. Gómez-Limón et al. investigated the priority of agricultural water resources allocation in different watersheds under drought conditions, and indicated that the priority can effectively improve water use efficiency [10]. In addition, similar studies were carried out on the administrative unit and irrigation district scale. For example, Li et at. determined the water allocation priority among different cities considering different socioeconomic development levels and national or regional policies [11]. Zhang et al. evaluated the priorities of 17 irrigation districts based on fairness, production efficiency and economic benefits, and then input them into a water resource optimization model to obtain a water resource allocation strategy [12]. However, the priority assessment system contains many spatial variability indicators, such as population, crop yield, etc. In addition, due to the limitations of data collection, it is difficult for current studies to consider smaller scales than administrative units and irrigation areas, which ignores the spatial variability of the system.

With the development of “3S” technology (i.e., remote sensing, Global Position System, and geographic information systems), more spatial information can be collected [13]. Some remote sensing data such as precipitation, evapotranspiration, and land use can be directly used in water resources management practices [14,15]. For example, remote sensing data have been utilized to evaluate the efficiency of water resources utilization, and to help optimize water resources allocation [14,16]. Moreover, some socioeconomic statistics can be estimated by appropriate methods and remote sensing data, which are helpful for water resources management [17]. The spatial unit of the collection boundary for socio-economic statistics is usually the administrative unit boundary. There are several limitations in spatial analysis and interdisciplinary application, including low spatial resolution and a mismatch between administrative regions and geographical units. Compared with traditional socio-economic data, the spatialization of statistical data has the following advantages: (1) Reflecting the spatial distribution characteristics of socio-economic data within the statistical region; (2) facilitating spatial analysis through grid-based socio-economic data; and (3) breaking the limitations of administrative boundary. Liang et al. used night-time light and land use data to spatialize the GDP of a Chinese city via a linear regression model, which can provide a reference for future urban planning and social and economic development [18]. Zhao et al. utilized different regression models to spatialize the GDP of South China and analyzed the performance of different models [19]. Sutton used night-time light data to map urban population density distribution and analyzed the trend of urban population density in the United States [20]. However, most relevant studies focused to gather the spatialization of GDP and population statistical data, few attempts have been made on agricultural economic statistical data, such as grain yield and agricultural output value [21]. Moreover, there are few reports on the application of statistical data spatialization methods in agricultural water management. Therefore, this research attempts to construct a spatialization method of agricultural output value and grain yield. After using the obtained spatialized statistical data and remote sensing data, we can ensure that the priority evaluated on a smaller scale contains more spatial information.

To formulate specific planning schemes, the priority of water users was integrated with optimization models. Generally, there are conflicting allocation objectives among different users including economic, social, and environmental impacts [3]. Thus, multiobjective models were developed for tackling multilateral interests from decision-makers [22,23,24]. In practical water-allocation problems, economic benefits, crop yield, and environmental impacts were the most common objectives, while fairness in water allocation was attracting more and more attention [25,26]. Zhang et al. used the Gini coefficient in water distribution to quantify the fairness among different water use sectors [23]. Li et al. measured the fairness of water distribution by the minimum difference of water use per unit area [3]. These indicators were presented as a function describing the absolute individual differences, such that penalties were incurred for any differences in individual distribution results (whether a customer is worse off or better off) [27]. Through such methods, interactions between allocation individuals, such as jealousy and comparison, can be hardly considered and reflected, which may reduce the resource allocation efficiency [28]. Envy is a measure that considers the differences in service quality between all possible pairs of customers [29]. This fairness index, based on the jealousy value, not only considers the psychology of comparison among the allocation objects, but also effectively avoids the loss of allocation efficiency. Therefore, this indicator is applied to the allocation of limited resources to consider the envy among the allocation objects in the allocation process.

In order to effectively manage agricultural water resources, an evaluation‒optimization framework has been developed. First, a water allocation priority evaluation system is constructed, and remote sensing data is involved. After applying this evaluation model at the canal scale, both the accuracy and spatial variability can be guaranteed in the evaluation results. Then a nonlinear multi-objective optimal allocation model of agricultural land and water resources (NMOLW), considering allocation efficiency and equity, was established. The presented method improves the existing statistical data spatialization method and is applied to the constructed water allocation priority evaluation system. Then we input the priority evaluation results into the NMOLW model, and obtained the agricultural water and land resource allocation of irrigation area plan by solving the model. This can: (1) Obtain the water allocation priority through considerations of the social, economic, and environmental factors at the canal scale; (2) improve the spatialization method of agricultural output value (specifically referring to the income from the planting industry) and grain yield statistical data by introduce normalized difference vegetation index (NDVI) and remote sensing data on potential yield. Potential yield of farmland refers to the farmland production potential estimated by the GAEZ (Global Agro Ecological Zones) model, considering soil, weather, the environment, and other factors [30]; and (3) illustrate the tradeoff between conflicting objectives, such as economic benefits, water supply priority, grain yield, and minimum envy. The method is applied to a practical case of an irrigation district in Northwest China and the optimization results can provide decision-making support for managers.

2. Development of Methodology

2.1. Overview of the Problem

In agricultural production, under the constraints of limited agricultural land and water resources, formulating more effective resource utilization strategies requires us to determine the water supply priority of competing water sectors. However, the current research on the priority of water supply cannot meet the requirements of precise management of water supply in practice. For example, most previous evaluation programs were carried out on administrative units and irrigation districts, which could not take into account the actual distribution of the canal system; the previous evaluation system does not consider the spatiality of the evaluation index, which has also limited the accuracy of the evaluation plan. Moreover, when allocating limited agricultural water resources, blindly pursuing distribution efficiency, such as economic benefits and grain output, will lead to an unfair distribution of water resources. Thus, to use land and water resources in a more efficient and sustainable way, decision makers are faced with the following issues: (1) How to evaluate the priority of water distribution on a smaller scale according to actual needs? (2) How to resolve the conflict between the boundary of statistical data and the boundary of the new evaluation unit? (3) How to balance distribution efficiency and distribution fairness in water resource distribution? Especially in the case of water scarcity, managers should not only consider food production and economic benefits, but also consider the fairness of distribution and user satisfaction with the water supply.

Therefore, the framework developed in this study (Figure 1) attempts to solve the above problems, and can be divided into three parts: (1) Spatialize the statistics of agricultural output value and grain yield; (2) construct a water allocation priority evaluation system on the canal scale; and (3) formulate a NMLOW model to obtain the irrigation water allocation.

2.2. Water Allocation Priority Evaluation

2.2.1. Spatialization of Statistical Data

Remote Sensing Data Preparation

The potential yield of farmland and land use (based on Landsat 8 remote sensing image, generated by human visual interpretation) remote sensing data were obtained from in the Resource and Environmental Science and Data Center of the Chinese Academy of Sciences (http://www.resdc.cn/, accessed on 20 May 2020) at a 1 km resolution. We used projection transformation tool to process the potential yield of farmland and land use data, and removed invalid values to obtain these data with a spatial resolution of 1 km in the WGS-1984 coordinate system. This study only focuses on agricultural output value (income from planting industry) and grain yield (yield of food crops). The fact that there will be no grain yield and agricultural output value without cultivated land needs to be stressed [17]. Therefore, the extracting tool in ArcGIS 10.2 (Environmental Systems Research Institute Inc., Redlands, CA, USA) was used to extract cultivated land, and we used the same tool to generate the potential yield of cultivated land.

Global 10 Days Synthesis of SPOT VEGETATION Images (VGT-S10) product was accessed via the Flemish Institute for Technology Research (VITO, https://www.vito-eodata.be, accessed on 11 June 2020) at a 1 km resolution. The time resolution of this product is 10 days, and we obtained 36 images in 2017. First, we did a projection transformation with the projection conversion tool. Then, the invalid value (negative value caused by the influence of clouds, water, ice and snow) was removed to obtain 36 periods of NDVI data within one year. Annual maximum NDVI generally represents NDVI at the height of the growing season, and reflects the growth of crops [31]. Through synthesizing 36 images into the annual maximum NDVI map with the maximum value composites (MVC) method, an annual maximum NDVI dataset for cultivated land was generated.

Regression Models

The multiple linear regression model is a very popular statistical regression model, that has been widely used in the spatialization of GDP, grain yield and other remote sensing applications [18,21,32,33]. However, in the spatialization of socio-economic data, a single regression model usually cannot accurately describe administrative units [19]. Therefore, more regression models such as quadratic polynomial models, power function models and exponential models have been introduced into remote sensing applications [34,35,36]. To examine the applicability of these models, we used them for measuring the quantitative relationships among the annual maximum NDVI, potential yield of farmland, statistical data of grain yield, and agricultural output value at the city level. The regression model performance was evaluated by goodness of fit value (R²) and the significance test value F [19,36,37].

Modification Model

The above models were established with the statistical data of city-level units, so when estimating the statistical data of each grid, there will be a large deviation [19]. To solve this problem, it is necessary to make amendments to the simulated data for each grid in the administrative unit [19,38,39]. Therefore, taking into account regional differences, city-level statistics are used for calculating the correction coefficient. The modified model is shown as follows:

$$A{P}_p=A{P}_i\times \frac{A{P}_t}{A{P}_s}$$

(1)

$$Y_p=Y_i\times \frac{Y_t}{Y_s}$$

(2)

A list of the nomenclatures for variables and parameters is provided in

Table 1. Nomenclature for variables and parameters.

Parameters	Definition
i	Index of response unit (i = 1, 2, 3…I)
c	Index of crop type (c = 1, 2, 3…C)
t	Index of month, from March to October (t = 1, 2, 3…T)
APp	The revised grid agricultural output value, CNY
APr	The grid agricultural output value estimated by regression model, CNY
APt	The agricultural output value counted by city units, CNY
APs	The sum of grid agricultural output value estimated by regression model within the scope of city unit, CNY
GYp	The revised grid grain yield, kg
GYr	The grid grain yield estimated by regression model, kg
GYt	The grain yield counted by administrative units, kg
GYs	The sum of grid grain yield estimated by regression model within the scope of each city unit, kg
Og	Measured value of the county g (provinces govern cities and cities govern counties)
Sg	Simulation value of the county g
ET0	Reference evapotranspiration, mm
Δ	Slope of the saturated vapor pressure and air temperature curve, kPa/°C
Rn	Net radiation, MJ/(m2·d)
G	Soil heat flux, MJ/(m2·d)
γ	Hygrometer constant, kPa/°C
T	Average air temperature, °C
u2	Wind speed at a height of 2 m above the ground, m/s
es	Air saturated vapor pressure, kPa
ea	Actual vapor pressure, kPa
Kc	Crop coefficient
ETict	Actual evapotranspiration of crop c in the t month of response unit i, mm
ETmct	Maximum evapotranspiration of crop c in the t month, mm;
λict	Water sensitivity parameters of crop c in the t month of response unit i
Aic	Irrigated area of crop c in response unit i, ha
Bc	The unit price of crop c, CNY
GPmic	Maximum yield per unit area of crop c in response unit i, kg/ha
Mai	Actual water consumption of response unit i, m3
Maei	The amount of water required by response unit i, m3
Mak	Actual water consumption of response unit k, m3
Mek	The amount of water required by response unit k, m3
hict	Available water in root zone of crop c in the month t of response unit i, mm
mict	Irrigation water of crop c in the month t of response unit i
pt	Effective precipitation of month t, mm
Geict	Groundwater recharge of crop c in the month t of response unit i
hm	Maximum available water in root zone, mm
Qt	Available water in month t, m3
η	Irrigation water utilization coefficient
Amaxic, Aminic	Maximum and minimum irrigation area, respectively, of crop c in response unit i, ha.
ZYminc	Minimum yield of crop c, kg
Ge	Total groundwater recharge, m3
p	Total effective precipitation, m3
α, σ	Proportional fairness coefficient, the degree of proportional fairness
x1, x2	Potential yield of farmland, kg/ha, NDVI

, which will be true for all equations.

Model Validation

The grain yield and agricultural output value of each administrative unit are calculated through ArcGIS based on the above model. The error between the simulation results and the statistical data is analyzed using county-level statistical data, and the performance of the model can be tested. The model’s accuracy was tested by mean absolute percent error (MAPE) and root mean square error (RMSE) [18]:

$$MAPE=\frac{1}{G}{\displaystyle \sum _g^G|\frac{O_g-S_g}{O_g}}|$$

(3)

$$RMSE=\sqrt{\frac{1}{G}{\displaystyle \sum _g^G{\left(O_g-S_g\right)}^2}}$$

(4)

2.2.2. Evaluation Index

To pursue fair and efficient water allocation among competing users, it is necessary to give attention to systems thinking and the integration of hydrologic, environmental, societal, and economic considerations [9]. Nevertheless, for agricultural water allocation, this issue is more specific, because agricultural water resources are mainly used for agricultural production. With the aim of being conducive to agricultural production, based on previous studies [9,10,12], the evaluation indices can be classified into three types: (1) Agricultural production conditions, such as land conditions, infrastructure construction, and regional policies; (2) economic benefits, such as agricultural output value and grain yield; and (3) equity, such as population and regional development level.

2.2.3. Evaluation Method

Division of Water Supply Response Unit

Canal systems are vital facilities in irrigation water distribution and, thus, it is more reasonable to divide water supply response units according to the canal system. Existing applications show that agricultural water management at the canal system scale can effectively improve management accuracy [40]. Thus, this research attempts to divide the water supply response unit by considering the distribution of the canal system and the current situation of land use and management.

The AHP-TOPSIS Evaluation Approach

As shown in Figure 2, the analytic hierarchy process (AHP) has a hierarchical structure that adapts to decision-making problems, so it has been widely used in multi-criteria decision-making. Using the experience and knowledge of local managers and AHP method can allow us make qualitative judgments on different objectives, and the specific method can be found in [41]. TOPSIS is a method to determine the closest positive ideal solution and the farthest negative ideal solution in multidimensional computing space [42]. Combining the AHP and TOPSIS methods can effectively measure the qualitative and quantitative problems in the evaluation and help us to obtain the ranking of different programs. The specific steps of the AHP-TOPSIS method can be found in the Appendix A.

2.3. Nonlinear Multi-Objective Optimization Model of Agricultural Land and Water

2.3.1. Crop Water Requirements

Crop water requirements are the main reference when formulating optimization models to guide water allocation [14]. The FAO56 PM equation recommended by the Food and Agriculture Organization of the United Nations (FAO) is the most commonly used method of traditional ground monitoring of ET₀ and Kc, and the maximum evapotranspiration of crops can be calculated according to the crop coefficient method recommended by the FAO [43]:

$$E{T}_0=\frac{0.408 \Delta {(\mathrm{R}}_n-\mathrm{G})+\gamma \frac{900}{T+273}{u}_2{(\mathrm{e}}_s-{\mathrm{e}}_a)}{\Delta +\gamma (1+0.34u_2)}$$

(5)

$$E{T}_m=K_c\times E{T}_0$$

(6)

2.3.2. Objectives

(1) Economic benefit

The Jensen model is a multiplicative crop water production function that can predict crop yield, and reflect the impact of water shortage on crop yield at various growth stages [44]. Therefore, this study uses this model to determine crop irrigation schemes.

$$Max={{\displaystyle \sum _i^I{\displaystyle \sum _c^C{\displaystyle \sum _t^T{\displaystyle \prod (\frac{E{T}_{ict}}{E{T}_{mct}})}}}}}^{{\lambda }_{ict}}\times A_{ic}\times GP{m}_{ic}\times B_c$$

(7)

(2) Grain yield

Irrigation districts often undertake the task of grain production. While pursuing economic benefits, grain yield should also be guaranteed.

$$Max={{\displaystyle \sum _i^I{\displaystyle \sum _c^{C-1}{\displaystyle \sum _t^T{\displaystyle \prod (\frac{E{T}_{ict}}{E{T}_{mct}})}}}}}^{{\lambda }_{ict}}\times A_{ic}\times GP{m}_{ic}$$

(8)

(3) Water allocation priority satisfaction

When water resources are limited, managers need to decide the water supply priorities of different water-using departments and direct the water supply under the guidance of this priority.

$$Max=W_i\times \frac{M_{ai}}{M_{ei}}$$

(9)

(4) Equity

In agricultural water allocation, fairness is seldom considered from the perspective of farmers’ psychological affect. Envy is a feeling of wanting what others have. In the model, it is defined as a user feeling unsatisfied compared with other users. When other users are less dissatisfied r, jealousy occurs [27]. Based on Sunarin’s research into the envy fairness index, this paper introduces the minimum envy fairness [27]. Envy value specifically refers to the difference between the unsatisfied water supply of the response unit i and that of the response unit k, namely:

$$Min={\displaystyle \sum _i^I{\displaystyle \sum _k^K\max }\{0,}(M_{ei}-M_{ai})-(M_{ek}-M_{ak})\}$$

(10)

2.3.3. Constraints

(1) Soil water balance constraint

The crop evapotranspiration is directly determined by the soil water content in the root zone. The soil water balance in the root zone can be simplified as follows [45]:

$$h_{ict+1}=h_{ict}+m_{ict}+p_t+G{e}_{ict}-E{T}_{ict}$$

(11)

$$h_{ict}\le h_m$$

(12)

(2) Evapotranspiration constraint

In practice, crop evapotranspiration cannot exceed the maximum evapotranspiration, nor can it be lower than the minimum evapotranspiration that satisfies crop growth [15].

$$0.6E{T}_{mct}\le E{T}_{ict}\le E{T}_{mct}$$

(13)

(3) Available water constraint

Irrigation water consumption cannot exceed the water available in the area.

$${\displaystyle \sum _i^I{\displaystyle \sum _c^C{A}_{ict}\times m_{ict}\le }}{Q}_t\times \eta$$

(14)

(4) Irrigation area constraint

The irrigated area cannot exceed the planting area, nor can it fall below the lower limit given by the manager.

$$A_{\min ic}\le A_{ic}\le A_{\max ic}$$

(15)

(5) Crop yield constraint

In order to meet the development requirements of the irrigation area, managers often set the output requirements of different crops in the irrigation area based on the production situation in previous years.

$${{\displaystyle \sum _i^I{\displaystyle \sum _t^T{\displaystyle \prod (\frac{E{T}_{ict}}{E{T}_{mct}})}}}}^{{\lambda }_{ict}}\times A_{ic}\times GP{m}_{ic}\ge Z{Y}_{\min c}$$

(16)

(6) Proportional equity constraint

In resource allocation, if the pursuit of minimum envy fairness is one-side, it will inevitably lead to neglect of users with less demand. To solve it, this study learns from the idea of proportional equity [46], giving each user the minimum requirement satisfaction.

$$\frac{M_{ei}}{M_{ai}}\ge \alpha$$

(17)

$$M_{ei}={\displaystyle \sum _c^C{\displaystyle \sum _t^{T}E{T}_{mct}}}\times A_{ic}$$

(18)

$$M_{ai}={\displaystyle \sum _c^C{\displaystyle \sum _t^{T}E{T}_{ict}}}\times A_{ic}$$

(19)

$$\alpha =\sigma \times \frac{Q+P+Ge}{M_e}$$

(20)

2.3.4. Solution Method

The weighted minimum deviation method is used for transforming the multi-objective model into a single objective model, and then a solution can be generated. The essence of this method is to transform multi-objective programming into single objective programming by normalizing each objective and eliminating the scale impact of an objective function value [22]. The specific calculation method can be found in [22,47], and managers’ preference for different objectives are obtained by AHP.

3. Application

3.1. Study Area

Hetao Irrigation District (HID) is located in Bayannur City, Inner Mongolia Autonomous Region, China (40°19′–41°18′N, 106°20′–109°19′E, altitude 1039 m), and its geolocation is shown in Figure 3. HID is an important commodity grain production base in China. Corn, sunflower, and wheat are the main crops in HID, accounting for 85% of the total sown area. There is little precipitation with an average annual precipitation of 130–220 mm and high annual evaporation of more than 2000 mm [48]. The irrigation and drainage canal system in HID is well constructed, with 13 main canals and 12 main ditches, and the total length of the canal system is more than 16,800 km. This study focuses on 12 canals, which are the most important canals in the irrigation area and are located in the core of the irrigation area. Approximately 90% of the irrigation water of the HID comes from the Yellow River. According to the statistics, the average annual irrigation water from the Yellow River is about 4.84 × 10⁹ m³. However, due to the decrease in annual runoff of the Yellow River, the water distribution of HID is reduced, and the gap between supply and demand of irrigation water is widened [49]. Water resources have become the main factor restricting the development of irrigation area. The challenges of the irrigation area are (1) how to manage water use on a canal system scale; (2) how to determine a reasonable priority for water distribution; and (3) how to allocate water resources in consideration of allocation efficiency and fairness.

3.2. Evaluation Index Selection

Based on Section 2.2.2 and combined with the actual situation of the irrigation area and data collection, seven evaluation indexes were selected in this study: Main canal length, canal design flow, number of branch canals, population, sowing area, agricultural output value, and grain yield.

3.3. Data

This study selected 2017 as the current year. The population, potential yield of farmland, and land use data from 2017 (spatial resolution: 1 km) were downloaded from the resource and environment data center of the Chinese Academy of Sciences (http://www.resdc.cn, accessed on 20 May 2020), and are shown in Figure 4. The synthesized 2017 maximum NDVI is shown in Figure 4b. The statistical data of agricultural output value and grain output of 12 cities in Inner Mongolia in 2017 are from the Inner Mongolia Statistical Yearbook of 2017 (Figure 5). The study area can be divided into 12 regions according to the main canal conditions: Yigan (YG, i = 1),Wulahe (WLH, i = 2), Yangjiahe (YJH, i = 3), Huangji (HJ, i = 4), Yongji (YJ, i = 5), Fengji (FJ, i = 6), Zaohe (ZH, i = 7), Shahe (SH, i = 8), Yihe (YH, i = 9), Tongji (TJ, i = 10), Changji (CJ, i = 11), Tabu (TB, i = 12). The meteorological data required by the FAO56 PM equation are all from the China Meteorological science data sharing service network (http://data.cma.cn, accessed on 10 May 2020). The canal system data of the 12 main canals (

Table 2. Canal system data of Hetao Irrigation District.

Irrigation Districts	Length of Main Canal (km)	Design Discharge of Canal (m3/s)	Number of Branch Canals
YG	17	65	17
WLH	53.4	27.4	15
YJH	58.5	53	25
HJ	81.4	58	34
YJ	49.4	106	39
FJ	105	60	36
ZH	52	29.2	23
SH	80	31.4	20
YH	83	33	24
TJ	67.9	31	24
CJ	53.5	24	20
TB	44.5	22	21

), the main crop planting structure (Figure 6), effective precipitation (

Table 3. Meteorological and hydrological data of Hetao Irrigation District.

Month	ETm (mm)			Effective Precipitation (mm)	Available Water (108 m3)
	Wheat	Corn	Sunflower
3	35.99	0	0	0	0
4	103.55	0	0	0	2.16
5	199.95	114.70	0	0	10.67
6	112.00	147.20	155.20	7.65	4.17
7	0	179.30	151.59	4.59	5.37
8	0	90.35	68.11	17.13	2.02
9	0	0	48.48	0	1.46

), and available water (Table 3) are all from the Hetao Irrigation District Management Bureau. The crop ETm data calculated based on the meteorological data are shown in Table 3. The crop yield and price data (

Table 4. Local price and minimum yield of crop.

Crop	Local Crop Prices (CNY/kg)	Crop Minimum Yield (104 kg)
Wheat	3.00	28,862.22
Corn	1.70	98,154.24
Sunflower	6.30	42,126.28

) are from the Inner Mongolia Statistical Yearbook of 2017. The statistical data of agricultural output value and grain yield of the five counties used for model accuracy verification are from the Bayannur City Statistical Yearbook of 2017.

4. Results and Discussion

4.1. Spatial Results of Statistical Data

4.1.1. Statistical Regression Model

When establishing the regression model of agricultural output value and grain yield, we used 12 cities in Inner Mongolia as the research objects (Figure 5). Then the correlation between remote sensing impact factors and statistical data can be analyzed. The results of the correlation analysis are shown in

Table 5. Correlation analysis results between remote sensing impact factors and statistical data.

	Impact Factors	Cultivated Land Area	NDVI	Potential Yield of Farmland
Statistical Data
Agricultural output value		0.93 **	0.78 **	0.65 *
Grain yield		0.87 **	0.86 **	0.73 **

Note: ** p < 0.05; * p < 0.1.

.

As shown in Table 5, there are a strong correlation among the three factors and agricultural output value, grain yield. All of them have passed the significance test. Among them, the cultivated land area has a best correlation with agricultural output value and grain yield. This reflects the principle that there is no agricultural output value and grain output without arable land is reasonable. In addition, the NDVI and potential yield of farmland have a good correlation with agricultural output value and grain yield, so these two factors can be used for predicting the statistical data. Existing studies have shown that, in the application of remote sensing, both the single-factor nonlinear model and the multi-factor multiple linear regression model have been applied well and achieved good estimation accuracy [18,19,21]. This study not only established a single factor model with NDVI and potential yield of farmland as variables, but also a multiple regression model including these two variables. Then a more suitable regression model can be selected by comparing the performance.

For comparison, the fitting results of all algorithms are shown in

Table 6. The fitting performance of different algorithms.

Factor	Algorithm	Agricultural Output Value		Grain Yield
		R2	F	R2	F
NDVI	Linear algorithm	0.61	15.37	0.74	28.203
	Quadratic polynomial algorithm	0.61	6.92	0.8	18.56
	Exponential algorithm	0.65	18.23	0.73	27.10
	Power function algorithm	0.72	25.42	0.77	33.124
Potential yield of farmland	Linear algorithm	0.42	7.223	0.53	11.18
	Quadratic polynomial algorithm	0.44	3.619	0.73	12.22
	Exponential algorithm	0.56	12.64	0.61	15.41
	Power function algorithm	0.63	16.87	0.56	12.56
Multiple linear regression		0.81	20.86	0.74	13.94

. Statistical regression results of agricultural output value show that the model established by NDVI is better than the potential yield of farmland. Among the four regression models, the power function model established by NDVI has the best fitting performance (R² = 0.7177, F = 25.42). Statistical regression results of grain yield also show that the model established by NDVI is better. Among these models, the quadratic function model and the power function model established by NDVI have better fitting results. In real cases, the agricultural output value and grain yield are not supposed to be negative. However, the linear and quadratic function models are prone to negative values, while the exponential function and power function algorithm would not be negative. Therefore, although the established linear model between NDVI and grain yield has a good fitting effect, we cannot select it due to the existence of negative values. In addition, considering the combined effects of NDVI and potential yield of farmland a multiple linear regression model is proposed. For preventing negative values, this study set the intercept of regression model to zero.

Then, a power function model of NDVI and agricultural output value, a power function model of NDVI and grain yield, a quadratic function model of NDVI and grain yield, and a multiple linear regression model were selected according to the comparison.

4.1.2. Spatialization Results

Five county-level administrative regions (Dengkou, Hangjin Houqi, Linhe, Wuyuan, and Urat Qianqi) in the study area were used as verification objects (Figure 3). The relationship between the spatial statistical results and the actual statistical results was analyzed to verify the accuracy. The results are shown in

Table 7. Simulation performance of regression models.

Factor	Algorithm	MAPE	RMSE
Agricultural output value (104 CNY)	Area weighted algorithm	28.37%	60,183
	NDVI power function algorithm	12.18%	40,283
	Multiple linear regression	25.54%	48,252
Grain yield (kg)	NDVI power function algorithm	16.69%	62,953
	NDVI quadratic polynomial algorithm	14.92%	61,161
	Multiple linear regression	14.66%	60,101

.

In order to demonstrate the superiority of the new model, the traditional area weighted algorithm is also used for agricultural output value spatialization (Table 7). Among the spatial models of agricultural output value, the precision of multiple linear regression and NDVI power function model are better than the area weighted algorithm. The NDVI power function model has the highest accuracy, and compared with the traditional method, the accuracy of this method is improved by 16%. This demonstrates that introducing impact factors can improve the simulation accuracy of the spatial model compared with the traditional algorithm. The error of the agricultural output value spatialization result obtained by the multiple linear regression algorithm is larger than that of the previous research [17]. However, previous studies included constants in the area where multiple regression models were established, which is contrary to the understanding that there is no agricultural output value without cultivated land. As for the spatialization of grain yield, previous studies established relationships between population density, land, and grain yield [50]. However, in fact, the correlation between population density and food production is not strong, which restricts the improvement of spatialization accuracy. Among the spatial models of grain yield, the multiple linear regression algorithm has the highest accuracy. According to the simulation results of grain yield, an algorithm having good fitting results does not necessarily mean high simulation accuracy. This also illustrates the necessity of choosing different regression models instead of a single regression model.

Based on the analysis, the NDVI power function algorithm was chosen for agricultural output value spatialization (Equation (21)), and the multiple linear regression was chosen for grain yield spatialization (Equation (22)). Through calculations, the statistical spatialization data was obtained. As shown in Figure 7, areas with higher agricultural output value are concentrated in the central region of the HID, including Hangjin Houqi, Linhe, and Wuyuan (Figure 7a). The areas with higher grain yield are Hangjin Houqi, and Linhe (Figure 7b). This information is consistent with the actual data, which also shows the validity of the simulation results.

$$y=645.17x_1^{4.0692}$$

(21)

$$y=0.112x_1+59.307x_2$$

(22)

4.2. Water Allocation Priority

According to the canal system distribution and control area, the whole irrigation area can be divided into 12 water supply response units. The division results are shown in Figure 8a. The water allocation priority of water supply response unit is evaluated by the AHP-TOPSIS method. The ranking results are displayed in Figure 8b. The darker color denotes the higher priority of water supply. The priority of the central irrigation area is significantly higher (especially in YJ, HJ and FJ response units) than others, reminding managers to pay more attention to the interests of the central area of HID. The agricultural output value and grain yield can be found in Section 4.1. Finally, the proximity of each response unit is calculated and the water allocation priority was obtained. The calculation results of proximity are shown in Figure 9.

For verifying the rationality of this priority evaluation, the average water diversion values of the main canals are compared with the evaluation results. The top five, YJ, HJ, FJ, YH, and YJH are selected, and the average water diversion volume of them in recent years is calculated. The statistical result of water amounts shows that YJ > HJ > FJ > YJH > YH, which is roughly consistent with the evaluation results. However, the situation in YJH and YH is somewhat different. This may be due to YJH’s larger grain planting area, and its water demand being higher than that of YH. From Figure 9, it can be noted that the proximity between YJH and YH is extremely close (YH = 0.522, YJH = 0.519). According to the actual situation, there is little difference between the two response units. On the basis of the evaluation and analysis, this study assumes that YJH and YH should have the same water allocation rights.

4.3. Optimization Results Analysis

The weight of each objective was obtained by AHP, and the results are shown in

Table 8. Weight of each objective.

Objective	Weight	CR
Economic benefit	0.18	0.04
Water allocation priority satisfaction	0.24
Grain yield	0.37
Equity	0.21

. When applying this method, the consistency ratio is 0.04 (less than 0.1), which indicates that the judgment matrix is a consistency matrix. The necessary data are input into NMOLW. Using the model solution in Section 2.3.4, the optimal schemes of water and land resources in HID were obtained.

Figure 10 shows the difference between irrigated area and planting area, indicating that corn planting areas can be well irrigated (Figure 10b) while some planting areas of wheat and sunflower are not able to be irrigated (Figure 10a,c, respectively), and the irrigation water for most sunflower area cannot be ensured. Existing crop irrigation water allocation schemes indicate that corn always gets more irrigation water under different available water conditions [51]. Corn has a greater planting advantage in HID, because corn yield and economic benefit are higher than those of wheat. Sunflowers have more non-irrigated areas than wheat, because managers pay more attention to food production. However, the reduction of guaranteed sunflower irrigation area will decrease farmers’ economic income. Therefore, when the irrigation area requires to ensure the grain yield under water shortage, it is needed to compensate farmers’ economic losses. Moreover, other water sources such as groundwater can be used temporarily for irrigation. Moreover, it can be found that HID currently has a large sunflower planting ratio, but the existing water resources cannot guarantee irrigation in so many areas. Therefore, in order to meet the needs of food production, this study suggests that HID should appropriately reduce the proportion of sunflower planting.

Figure 11 shows the changes in irrigation quota before and after optimization, indicating that the irrigation quotas for corn and sunflower have increased while the irrigation quotas for wheat have decreased. This proves that wheat planting is indeed inferior to corn and sunflower. In all regions, WLH’s and YJH’s irrigation quotas have the highest level of improvement. In fact, these two areas are located in very important areas in the HID [40], and the current low irrigation quota cannot give full play to the local production potential. Unlike in other regions, the irrigation quota of all crops in FJ has been reduced, which shows that the irrigation quota in this area can be reduced under the premise of ensuring agricultural production. In addition, YJ and HJ with the highest priority have the lowest increase in irrigation quota levels. This shows that these two areas have received attention already, so it is enough to maintain the current level of irrigation.

When the weighted minimum deviation method is applied to solve the multi-objective model, the final weighted coordination satisfaction of each objective is 0.82 (the closer to 1, the better the coordination results between the objectives), indicating that tradeoffs are well addressed. Comparing the optimization results with the results from 2017, the grain yield and economic benefits after optimization have been improved by 19.8% and 4%, respectively. This indicates that the optimal allocation of agricultural water and land resources is helpful to agricultural production. In fact, HID has increased grain production capacity in recent years to complete grain production tasks, but the optimization results prove that there is still room for improvement in grain production capacity. Therefore, to reduce the pressure on food production, managers can further increase production capacity based on the optimization results.

4.4. Discussion

Although the introduction of NDVI and farmland production potential data improved the accuracy of agricultural output value and grain yield estimation, the results of spatialized data still contained some uncertainties, caused by the following factors. Firstly, the reliability of the statistical data directly affects the accuracy of regression models and spatialized results since the statistical data are the basis for the regression analysis [19]. Secondly, despite the preprocessing procedures used in NDVI and the potential yield of farmland data to reduce background noise, the corrected data may still contain some noise, which can affect the statistical data estimation. Thirdly, using NDVI and farmland production potential data only, it is difficult to fully and accurately express the spatial heterogeneity of agricultural output value and grain yield distribution, which may lead to overestimation or underestimation at the local area.

Table 9. Accuracy verification results of each county.

County	Grain Yield (103 kg)			Agricultural Output Value (104 CNY)
	Statistical Data	Estimation Results	MAPE (%)	Statistical Data	Estimation Results	MAPE (%)
Linhe	557,750.00	517,682.00	7.18	346,934.50	280,051.00	19.28
Wuyuan	419,200.00	494,521.00	17.96	244,794.50	298,946.00	22.12
Hangjing Houqi	420,180.00	417,076.00	0.74	248,715.20	225,567.00	9.31
Urat Qianqi	500,225.00	413,854.00	17.23	317,690.50	305,545.00	3.82
Dengkou	191,000.00	248,555.00	30.13	78,170.99	83,150.00	6.37

shows the accuracy verification results of the five counties. It can be seen that the MAPE in most areas is less than 20%, and the simulation accuracy can be accepted. However, the simulation results in Wuyuan and Dengkou have large errors, which may have been caused by the impact of the planting structure. For areas where the proportion of food crops is high, the agricultural output value will be overestimated (Wuyuan). For areas dominated by cash crops, the grain yield will be overestimated (Dengkou). Finally, many publications indicate that NDVI does have a close relationship with agricultural production. However, the existing relationships are based on statistical relationships, not on the formation mechanism of crop yields [21,50]. Gross primary productivity is highly correlated with biological productivity and commonly used in crop yield estimations [52]. Therefore, future studies should explore the use of gross primary productivity. Overall, this is likely to enhance the accuracy of the estimation by improving the availability and quality of statistical data and remote sensing data, introducing multi-source data, such as planting structure data.

To demonstrate the superiority of minimum envy fairness, this study used the Gini coefficient as a fairness index in optimization. Single-objective models with minimum envy fairness (Equation (10)) and Gini coefficient (Equation (23)) are established, which have the same constraints as the NMOLW model. Water supply satisfaction, economic benefit, and grain yield are selected as indicators for assessing the performance of the two models. The calculation function of water supply satisfaction is shown in Equation (24):

$$\min Gini=\frac{{\displaystyle \sum _i^I{\displaystyle \sum _k^K\left|\frac{M_{ai}}{A_i}-\frac{M_{ak}}{A_k}\right|}}}{2I{\displaystyle \sum _i^I\frac{M_{ai}}{A_i}}}$$

(23)

$$p={\displaystyle \sum _i^I\frac{M_{ai}}{M_{ei}}}$$

(24)

As shown in

Table 10. Performance of two equity indices.

Attribute	Gini	Minimum Envy Fairness
p	0.39	0.47
Grain yield (108 kg)	18.07	21.83
Economic benefit (108 CNY)	74.43	80.98

, minimum envy fairness results in higher food production and economic benefits. When the optimization results of Gini coefficient are used they show that the greater the water demand, the higher the degree of dissatisfaction. Moreover, this study suggests that the minimum envy fairness should not be used alone when a significant water demand difference exist among users, but should be combined with other methods for constructing the fairness index. Minimum envy fairness well reflects the comparison and jealousy among water users. In addition, the minimum envy fairness can also be used to consider the distribution from the perspective of demand, avoiding efficiency loss. The Gini coefficient cannot consider these factors in practical resources allocation problems [27].

With the help of spatialized statistical data, the evaluation of water supply priority can be carried out on different scales, leading to higher accuracy than the previous evaluation. Therefore, in future research, we can extend the spatialization of statistical data to more uses such as agricultural industry analysis and agricultural regional planning. Due to the limitation of selection factors, the spatialization method of grain yield constructed in this paper is more suitable for the main grain producing areas planted as one-year crops. Therefore, how to measure the impact of different cropping systems on the spatialization of agricultural economic statistics is also a topic that we need to give attention to in the future. In addition, the NMOLW model pays more attention to socio-economic goals, but involves insufficient consideration of ecological aspects. Moreover, because of the limitations of the data collection, the factors considered in the evaluation of water allocation priority are limited, which affects the evaluation accuracy. These should be addressed in future work.

5. Conclusions

To determine the water allocation priority and manage agricultural water effectively, this study carried out the water allocation priority evaluation, and then a NMOLW model was proposed for considering fairness and efficiency. The water supply response units were divided according to the practical water supply condition of the study area. Furthermore, in order to solve the boundary conflict between the response unit and the original statistical data, a spatial method of analyzing the statistical data was introduced. This study has the following advantages:(1) Improving the spatial simulation accuracy agricultural output value and grain yield; (2) measuring the precise priorities including more spatial information; (3) fully considering the tradeoffs between allocation efficiency and equity; and (4) providing decision-making support for irrigation districts to formulate water use plans.

To illustrate the applicability of the proposed approach, this method has been applied to HID. The results show that: (1) the spatialization method of grain yield and agricultural production value achieved good simulation accuracy in HID; (2) we can determine the priority of water allocation at the scale of the canal system in HID; and (3) there is still room for increasing grain output in HID. The approach proposed in this study is expected to support sustainable agricultural water management in similar regions.

In fact, due to the variability of water supply response unit division and evaluating indexes, the method can be changed to suit different regions. There are many factors affecting the spatial distribution of agricultural output value and grain yield statistics, such as the detailed planting structure. These spatialization methods can be further improved.