CN113177345B - 一种网格化作物种植布局优化方法 - Google Patents

Abstract

本发明提供了一种网格化作物种植布局优化方法，该方法通过对目标种植区域建立网格模型，获取每个网格单元内的作物潜在单产和适宜种植面积、作物最小灌溉需水量、作物施肥强度、耕地面积以及灌溉面积数据，结合目标种植区的县级统计种植数据，分配得到每个网格单元内的作物种植数据，以灌溉需水总量最小为优化目标，基于预定的约束条件进行种植布局优化，得到优化后各个网格单元内的作物种植面积。本发明的优化方法是基于网格单元进行，得到的优化布局结果更加精准详细，可以科学指导作物种植布局调整，有效缓解地下水位下降等生态环境问题。

Description

一种网格化作物种植布局优化方法

技术领域

本发明属于种植布局优化调整技术领域，具体涉及一种网格化作物种植布局优化方法。

背景技术

目前，基于水资源展开的作物种植布局优化研究很多。研究人员多通过约束地区水资源消耗最少来调整得到节水的最优布局。然而，现有的优化布局存在一定的缺陷。第一，当前的优化布局多在区域尺度上进行，而区域尺度上的研究具有较大的空间不确定性。也有少部分研究将作物布局精细到网格上，然而优化方法多为线性优化，优化效果有待提升。第二，当前的优化研究中多对作物产量进行限制，未考虑作物不同的减产情景。第三，当前多数研究通过改变作物种植制度，如两年三熟，三年四熟等来达到节水目标，而改变种植制度需要改变农户的种植意愿，实施起来较为困难。因此，通过设定不同的作物产量缩减目标，结合非线性优化方法，在考虑农民种植习惯和灌溉设施等实际情况下，对作物的种植布局开展网格尺度上的优化调整，能够更加有针对性的缓解重点地区地下水位下降等生态环境问题，同时保证调整措施可以有效实施。

发明内容

本发明所要解决的技术问题是提供一种网格化作物种植布局优化新方法。在模拟过程中，通过设定不同作物减产比例，采用通用代数建模系统(The General AlgebraicModeling System，GAMS)提供的Python API接口，选择Conopt4求解器中的非线性优化求解算法进行优化模拟。

根据本发明的一方面，提供了一种网格化作物种植布局优化方法，包括：步骤S1、对目标种植区域建立网格模型，获取每个网格单元内的作物潜在单产和适宜种植面积；步骤S2、获取每个网格单元内的作物最小灌溉需水量；步骤S3、根据目标种植区的县级尺度作物实际种植数据，分配得到每个网格单元内的作物种植数据，种植数据包括灌溉和雨养种植条件下的作物种植面积和单产数据；步骤S4、获取每个网格单元内的作物施肥强度数据；步骤S5、获取网格单元内的耕地面积和灌溉面积；步骤S6，以灌溉需水总量最小为优化目标，基于预定的约束条件和前述步骤得到的数据进行种植布局优化，得到优化后各个网格单元内的作物种植面积。

其中，所述步骤S2进一步包括：步骤S21，计算网格单元内的参考作物潜在蒸散量；步骤S22，结合参考作物潜在蒸散量、作物系数和水分胁迫系数，计算生育期内作物潜在蒸散量和作物实际蒸散量；步骤S23，计算网格单元的有效降雨量；步骤S24，基于作物潜在蒸散量和有效降雨量计算每个网格单元的作物最小灌溉需水量。

进一步，所述步骤S3包括：步骤S31，计算县级尺度灌溉和雨养种植条件下作物的单产和播种面积；步骤S32，计算作物在灌溉和雨养种植条件下的适宜种植面积；步骤S33，计算作物的市场可达性；步骤S34，计算作物先验分布参数；步骤S35，基于最小化交叉信息熵方法求解作物在网格单元内分配的种植比例；步骤S36，计算每个网格单元内作物的单产和种植面积。

优选的，在步骤S6之前，还包括步骤S6-1，对步骤S1-S5获得的网格数据转化为点状数据。

优选的，所述步骤S6中，所述优化目标设定为优化后种植区作物灌溉需水总量最小；其中，灌溉需水总量为步骤S2中得到的每个网格单元内最小灌溉需水量乘以对应的每个网格单元优化后的种植面积得到每个网格单元内的灌溉需水总量，再将种植区范围内所有网格单元上的灌溉需水总量相加，即得到优化后种植区作物灌溉需水总量。

优选的，所述优化目标进一步表征为下述目标函数：

IRD_j＝∑_io_ijl×CropArea_jl×Dirri_i×10,l＝irri (6-1)

其中，IRD_j表示全种植区作物j的灌溉需水总量；o_ijl表示优化后网格单元i内作物j在灌溉I或雨养R种植条件l下的种植比例；CropArea_jl表示作物j在全种植区实际的灌溉I或雨养R种植条件l下的总种植面积；Dirri_i表示网格单元i上的最小灌溉需水量；10为转换系数。

优选的，通过下述8个约束条件实现所述优化目标：

(1)在每个网格单元上灌溉和雨养种植条件下分配的作物种植面积加和小于等于该网格单元上的耕地面积；

(2)每个网格单元上分配的作物灌溉面积小于等于该网格单元上的灌溉面积；

(3)每个网格单元上分配的作物灌溉面积或雨养面积小于等于网格单元上灌溉或雨养种植条件下的作物适宜性种植面积；

(4)每个网格单元上分配的作物雨养面积小于等于原布局上的作物雨养面积；

(5)优化后作物的种植区在原作物种植网格单元范围内；

(6)优化后作物的总产大于等于不同减产比例与优化前种植区作物总产的乘积；

(7)优化后种植区肥料施用总量小于等于优化前的总施用量；

(8)优化后种植区作物灌溉需水总量小于等于优化前的作物灌溉需水总量。

如上所述，本发明提出了一种网格化作物种植布局优化方法，用于缓解当前地区用水紧缺以及地下水位下降等生态环境问题。该方法中考虑了不同作物减产限制条件，分别模拟了作物总产不变、减少情景下以节水为目标的可能布局。

该方法中的优化算法是基于网格单元进行的，通过非线性优化求解，以地区作物的灌溉需水总量最小为目标进行种植布局优化。同时，优化建模中设定了八个约束条件，分别保证在每个网格单元上灌溉和雨养种植条件下分配的作物种植面积加和小于等于该网格单元上的耕地面积、每个网格单元上分配的作物灌溉面积小于等于该网格单元上的灌溉面积、每个网格单元上分配的作物灌溉或雨养种植面积小于等于网格单元上灌溉或雨养种植条件下的作物适宜性种植面积、每个网格单元上分配的作物雨养面积小于等于原布局上的作物雨养面积、优化后作物的种植区保证在原作物种植网格单元范围内、优化后作物的总产大于等于不同减产比例与优化前种植区作物总产的乘积、优化后种植区肥料施用总量小于等于优化前的总施用量。

通过层层约束，能够保证得到更加合理的优化布局结果，并且能够在网格单元上精准详细的提出作物缩减目标，从而达到节水目的。此外，该方法考虑了农民种植习惯和灌溉设施等实际情况，能够在不改变农户意愿的情形下调整作物种植布局，实施起来较为可行。

本发明的上述技术方案具有如下有益的技术效果：

本发明是一种在网格尺度上对作物种植布局进行优化的方法，该发明进一步完善和发展了现有优化布局方法的可行性等问题。该方法考虑了不同减产情景下作物种植布局的最佳节水方案。此外，由于该优化方法是基于网格单元进行，得到的优化布局结果更加精准详细，可为该地区农业供给侧结构性改革中作物优化布局调整提供科学依据，有助于国家或地方政府更加方便的展开有针对性的、具体的规划措施，从网格单元进行管理、调整，从而有效的缓解地下水位下降等生态环境问题。

附图说明

图1显示了本发明的网格化作物种植布局优化方法的流程示意图。

图2显示了本发明的网格化作物种植布局优化方法的技术原理示意图。

图3是根据本发明的布局优化方法实施后灌溉小麦种植面积空间分布图。

图4是根据本发明的布局优化方法实施后雨养小麦种植面积空间分布图。

图5是根据本发明的布局优化方法实施后小麦总种植面积空间分布图。

图6是根据本发明的布局优化方法实施后小麦灌溉需水量空间分布图。

图中阴影部分表示每个网格单元上的数值大小，颜色越深，表示该网格单元上的数值越大。

具体实施方式

为使本发明的目的、技术方案和优点更加清楚明了，下面结合具体实施方式并参照附图，对本发明进一步详细说明。应该理解，这些描述只是示例性的，而并非要限制本发明的范围。此外，在以下说明中，省略了对公知结构和技术的描述，以避免不必要地混淆本发明的概念。

各类面积定义：

作物种植面积：作物实际种植占用的耕地面积。

作物播种面积：考虑作物不同种植季节、种植次数的种植面积。如在一年两熟区，一块地一年内种植两次同种作物，则播种面积为两次种植面积的加和。

作物适宜性种植面积：作物在灌溉或雨养种植条件下适合种植的区域面积，为理论情况。该面积综合考虑了气候、土壤和地形等要素。

当前作物种植布局：是经过多年农户的种植，依据实际情况所形成的(实际情况)，该布局形成过程中不仅受到气候、地形等条件的影响，也受到人为因素影响，如作物适宜种植区被建筑用地或林草地占用等情况。

图1显示了本发明的网格化作物种植布局优化方法的流程示意图。

图2显示了本发明的网格化作物种植布局优化方法的技术原理示意图。

如图2所示，本发明提供了一种网格化作物种植布局优化方法。本发明中种植布局优化方法的技术原理主要包含三个阶段，第一阶段是优化布局建模所需的基础数据准备过程(包括步骤S1-S5)，第二阶段是优化布局方法的建模，通过模型计算输出优化布局(包括步骤S6)，第三阶段是将输出的优化布局结果空间化(步骤S7)。

在优化布局建模处理过程中，通过设定不同作物减产比例，采用通用代数建模系统(The General Algebraic Modeling System，GAMS)提供的Python API接口，选择Conopt4求解器中的非线性优化求解算法进行优化模拟。具体的，在模型中的优化函数及目标函数采用了非线性函数进行模拟。这里，Conopt4求解器是由外界内嵌连接到GAMS系统中的求解器，用于求解大型的非线性规划问题。在本发明的优化求解过程中，在GAMS系统可以直接调用该求解器。

下面以黄淮海区小麦种植为例详细介绍各个步骤的处理过程，但本发明不限制于此，本发明的优化模型和方法也适用于其他类似农作物的种植布局优化，例如水稻、大豆、玉米等，根据实际情况调整参数即可。

步骤S1、对目标种植区域建立网格模型，获取每个网格单元内的小麦潜在单产和适宜种植面积。

综合各类气象数据、土壤数据、土地利用数据、地形数据等，通过全球农业生态区划(GAEZ)模型模拟，对目标种植区域建立网格模型，划分成若干个网格单元，基于模型输入网格尺度数据，进一步得到种植区网格尺度上小麦在灌溉和雨养种植条件下的潜在单产和适宜性种植面积数据。该适宜性种植面积表示网格单元内气候条件、土壤条件等均适合小麦生长的区域面积。GAEZ模型依据每个土地评价单元上的气候(温度、降雨、辐射)、土壤条件等，测算小麦潜在单产和适宜小麦种植的面积。该结果为网格尺度上的结果，即每个网格单元内有对应的小麦在灌溉和雨养种植条件下的潜在单产(单位：kg/ha)和适宜种植面积(单位：ha)。

GAEZ模型共包含7个模块，分别用于农业气象数据分析和农业气象指标的编制、计算各种作物农业气象评估和雨养潜在生物量/产量、计算农业气象条件约束下的产量、计算土壤和地形限制下的产量、计算作物的潜在产量、实际作物产量数据处理和计算主要粮食作物和纤维作物的产量差。GAEZ模型根据每个土地评价单元上的气候条件和土壤条件，给出不同作物在该土地单元上种植的适宜性等级，而后根据研究时间段内的平均气候条件，对适宜的作物采用逐级限制法来计算土地粮食生产潜力，即按光合生产潜力(仅光照限制)—光温生产潜力(光照和温度限制)—气候生产潜力(光照、温度和水分限制)—农业-气候限制条件下的生产潜力(光照、温度、水分和病虫害、霜冻等灾害限制)—作物生产潜力(光照、温度、水分和病虫害、霜冻等灾害限制和土壤条件限制)，最后输出潜在产量。

该步骤中气象数据包含气温、降雨、风速、湿度、日照时数、气压等，来源于中国气象数据网；土壤数据来源于中国科学院南京土壤所以及HWSD土壤数据库，主要包含土壤有机碳、pH值、土壤持水能力、土壤深度、土壤阳离子交换能力、土壤质地等；土地利用数据来源于中国科学院资源环境科学数据中心；地形数据来源于地理空间数据云平台，地形数据指地面高程数据，即海拔高度。

步骤S2、获取每个网格单元内小麦的最小灌溉需水量。

步骤S21，计算每个网格单元内的参考作物潜在蒸散量。

应用联合国粮食与农业组织(FAO)推荐的以能量平衡和水汽扩散理论为基础的Penman-Monteith(P-M)公式，计算参考作物潜在蒸散量ET₀。

式中，Δ表示为温度曲线和饱和水气压相关的斜率，KPa/℃；R_n表示地表净辐射量，MJ/m²·d；G表示土壤热通量，MJ/m²·d，一般取0；γ表示干湿度常数，KPa/℃；T表示空气平均温度，℃；U₂表示2m处风速，m/s；e_s表示饱和水气压，KPa；e_a表示实际观测水气压差，KPa。

步骤S22，计算生育期内小麦潜在蒸散量和小麦实际蒸散量。

本步骤结合参考作物潜在蒸散量ET₀、作物系数K_c,t和水分胁迫系数K_s,t，计算生育期内逐日的小麦潜在蒸散量PET_c,t和实际蒸散量ET_a,t。

PET_c,t＝K_c,tET_0,t (2-2)

ET_a,t＝K_c,tK_s,tET_0,t (2-3)

公式(2-2)为作物在第t天潜在蒸散量PET_c,t的计算，其中K_c,t为第t天的作物系数，ET_0,t为第t天的参考作物潜在蒸散；

公式(2-3)为作物在第t天实际蒸散量ET_a,t的计算，其中K_s,t为第t天的作物水分胁迫系数，具体计算公式如下：

式中，s_t为作物根部平均深度的土壤含水量；S_max为根区最大土壤有效含水量；P为作物从根区土壤中吸收的无水分胁迫的水分含量占根区最大土壤有效含水量S_max的比例；

S_t＝S_t-1+Δt(P_eff,t+I_t-ET_a,t-R_t) (2-5)

式中，S_t-1为前一时期的土壤含水量，Δt为1天，P_eff,t为第t天的有效降雨，I_t为第t天的灌溉用水(仅针对灌溉作物，雨养条件下为0)；ET_a,t为第t天的作物实际蒸散量；R_t为第t天的径流；

R_t＝(P_eff,t+I_t)×(S_t/S_max)^γ (2-6)

式中，γ为常数，在雨养地区为2，灌溉地区为3；

P＝P_std+0.04(5-PET_c,t) (2-7)

式中，P_std为常数，PET_c,t为作物在第t天的潜在蒸散量。

步骤S23，计算网格单元内的有效降雨量。

公式(2-8)为有效降雨量P_eff,t的计算，其中P_t表示第t天的实际降雨量，当第t天的平均温度低于0℃时，认为ΔSnow等于P_t，当第t-1天有积雪，且第t天的平均温度高于0℃，则认为ΔSnow等于平均温度乘以4mm/℃。

P_eff,t＝P_t-ΔSnow (2-8)

步骤S24，基于小麦潜在蒸散量和有效降雨量计算网格区域的小麦最小灌溉需水量。

该最小灌溉需水量表示小麦生长过程中能够满足其自身生长发育的理论最小需水量，即除去降雨补充的水分后，仍然需要通过灌溉方式补充的水资源量。该步骤得到的小麦最小灌溉需水量为网格尺度上的结果。最小灌溉需水量Dirri的计算公式如下：

其中，当生育期内小麦潜在蒸散量∑PET_c,t大于等于有效降雨量∑P_eff,t时，最小灌溉需水量Dirri等于两者的差值；当生育期内小麦潜在蒸散量∑PET_c,t小于有效降雨量∑P_eff,t时，则最小灌溉需水量Dirri为0。

步骤S3、根据目标种植区的县级尺度小麦实际种植数据，分配得到每个网格单元内小麦的种植数据，种植数据包括灌溉和雨养种植条件下小麦的种植面积和单产数据。

获取种植区范围内各县域小麦的实际种植情况，包括播种面积、单产等数据，通过作物产量空间分布模型(SPAM)和最小化交叉信息熵方法，计算模拟得到网格尺度上灌溉和雨养种植条件下小麦的种植面积和单产数据。具体模拟算法如下的子步骤S31-S36。

本发明的优选实施例中，以最小化交叉信息熵和作物产量空间分布模型(Spatialproduction allocation model,SPAM)为理论基础，结合黄淮海区小麦种植的实际情况进行修正，得到适用于黄淮海区的小麦空间分布模拟模型，进而模拟实现小麦种植面积和实际产量的空间分布，模型的主要计算流程如下。

步骤S31，计算灌溉和雨养种植条件下小麦的单产和播种面积。

依据现有研究结果得到黄淮海区各流域范围内小麦在灌溉和雨养种植条件下单产的比值参数∝_j(下表1)。黄淮海区小麦灌溉和雨养种植条件下种植比例Percent_jl(j表示不同作物类型，l表示不同种植条件，即灌溉条件I和雨养条件R)则依据“2010年全球空间分类作物统计数据v2.0”的数据结果进行计算，该数据集中包含全球2010年42种作物在灌溉和雨养种植条件下的种植面积，灌溉和雨养种植条件下的种植面积与总种植面积(灌溉加雨养)的比值即为Percent_jl。依据小麦的灌溉比例Percent_jI、雨养比例Percent_jR和比值参数∝_j可以计算灌溉条件和雨养种植条件下的单产。

Yield_jI＝∝_j×Yield_jR (3-2)

其中，Yield_j为小麦的县级统计单产，Yield_jR表示县级尺度雨养小麦单产，Yield_jI表示县级尺度灌溉小麦单产，通过乘以比值参数∝_j可以得到；

HarvestArea_jl＝Percent_jl×HA_j (3-3)

其中，HarvestArea_jl表示灌溉或雨养种植条件下的小麦播种面积，HA_j表示县级统计数据的小麦总播种面积。

表1黄淮海区不同流域灌溉和雨养种植条件下小麦主要种植参数

步骤S32，计算小麦在灌溉和雨养种植条件下的适宜种植面积

研究通过步骤S1的全球农业生态区划(Global Agro Ecological Zone,GAEZ)模型模拟得到小麦在灌溉和雨养种植条件下的适宜种植面积数据。

步骤S33，计算小麦的市场可达性

本研究采用人口密度计算市场可达性。当前研究认为人口密度和市场可达性之间存在较高的相关性，因为作物生产者也是消费者，同时更高的人口密度意味着更大的市场潜力和更低的运输成本。因此可以依据下列公式计算市场可达性Access_ic：

其中，n＝1～N，MaxPop_c为县域c内最大的人口密度，MinPop_c为县域c内最小的人口密度，K_c为县域c内最大人口密度和最小人口密度差均分为100份的间隔。

步骤S34，计算小麦先验分布参数

基于上述公式得到的结果，依据公式(3-6)和(3-7)求得小麦生产的成本参数Cost_ijl以及先验分布参数p_ijl，该参数可以代表作物j(小麦)在网格单元i中出现的可能性大小，越高的成本参数意味着该网格单元中的分配比例越大。

Cost_ijl＝Access_ij×PotYield_ijl×Percent_jl (3-6)

其中，Access_ij为网格单元i内作物j的市场可达性，PotYield_ijl为网格单元i内作物j(小麦)在灌溉或雨养种植条件l下的潜在单产，Percent_jl为作物j(小麦)在灌溉或雨养种植条件l下的种植比例；p_ijl为网格单元i内作物j在灌溉或雨养种植条件l下分配的先验分布参数。

步骤S34是通过一些参数，如市场可达性、适宜性面积等得到网格尺度上的小麦先验分配参数，即先得到大致每个网格上小麦种植面积可能的分配比例p_ijl；然后通过步骤S35，在各个约束条件的限制下，得到最精确的每个网格上小麦种植面积分配的比例s_ijl，最后结合小麦的总种植面积乘以每个网格单元上的分配比例s_ijl，得到每个网格上的种植面积。

步骤S35，基于最小化交叉信息熵求解小麦分布比例

交叉熵用于描述两个可能分布之间的差异，可以反映系统的稳定性，交叉熵越小，变量的不确定性越小，系统越稳定。对作物分布进行求解时，通过在满足耕地、适宜种植面积、灌溉面积等约束条件下，最小化先验分布参数p_ijl和作物分布s_ijl之间的交叉熵，可以使s_ijl尽可能接近真实的概率分布：

En_ijl＝∑_i∑_j∑_ls_ijllns_ijl-∑_i∑_j∑_ls_ijllnp_ijl

其中En_ijl为交叉熵，s_ijl表示网格单元上分配的小麦在灌溉或雨养种植条件下的种植比例。

在进行最小交叉熵的模拟时，需要同时满足以下约束条件：

∑_is_ijl＝1 (3-8)

∑_lCropArea_jl×s_ijl≤lucc_i (3-9)

CropArea_jl×s_ijl≤SuitArea_ijl (3-10)

∑_i∈c∑_lCropArea_jl×s_ijl＝StArea_jc (3-11)

∑_lCropArea_jl×s_ijl≤IrriArea_i,l＝irri (3-12)

其中，IrriArea_i为网格单元i上的灌溉面积；lucc_i为网格单元i上的耕地面积；CropArea_jl表示作物j(小麦)在灌溉或雨养种植条件l下的种植面积；SuitArea_ijl表示网格单元i上作物j(小麦)在灌溉或雨养种植条件l下的适宜性种植面积；StArea_jc表示各县域c的统计数据中作物j(小麦)实际播种面积与作物j(小麦)熟制的比例；公式(3-8)约束了作物j(小麦)在种植条件l下，各县域内所有网格单元i上分配的种植比例加和等于1；公式(3-9)约束了作物j(小麦)在网格单元i上，在灌溉和雨养种植条件l下分配的种植面积加和小于等于网格单元i上的耕地面积；公式(3-10)约束了作物j(小麦)在网格单元i上灌溉或雨养种植条件l下分配的种植面积小于等于作物适宜种植面积；公式(3-11)约束了各县域内作物j(小麦)在所有网格单元i上分配的种植面积等于统计数据中作物j(小麦)实际播种面积与作物j(小麦)熟制的比例；公式(3-12)约束作物j(小麦)在网格单元i上分配的灌溉面积小于等于该网格单元i上的灌溉面积；

其中CropArea_jl的计算公式如下：

其中，CropIntensity_jl为作物j(小麦)在灌溉或雨养种植条件l下的熟制。

本发明的模型中，目标函数是存在对数关系的非线性问题，且需要优化求解大量的s_ijl变量。因此，可以通过通用代数建模系统(The General Algebraic ModelingSystem)提供的Python API接口，选择非线性函数编程实现黄淮海区的优化求解。

步骤S36，计算每个网格单元内小麦的单产和种植面积。

依据求得的s_ijl，可以计算各网格单元i上分配的小麦种植面积Area_ijl：

Area_ijl＝s_ijl×CropArea_jl (3-14)

其次，依据小麦在灌溉和雨养种植条件下的潜在单产PotYield_ijl可以求得各县域内的平均潜在单产Yavg_jl，并依据公式(3-16)计算得到各网格单元i上小麦在不同种植条件下的实际单产Y_ijl。

上述步骤S31-S36分为四个阶段：

第一阶段即S31分配灌溉和雨养种植条件下的小麦单产和播种面积。由于从统计年鉴搜集的数据是不分灌溉和雨养种植条件的，而后续算法的应用过程中需要划分，因此首先需要将统计数据按照灌溉和雨养种植条件进行分类。

第二阶段主要用于计算小麦分布的先验参数，包含步骤S32-S34。其中，步骤S32-S33为数据准备过程，步骤S34小麦先验分布参数的计算过程。

第三阶段步骤S35，主要依据小麦先验分布参数计算每个网格单元上最接近于真实分布的概率s_ijl。

第四阶段步骤S36，即依据步骤S35求得的s_ijl计算每个网格单元上小麦的单产和种植面积，即步骤S3的最终计算结果。

本发明步骤S3的算法原理即通过结合县域尺度上的统计数据，按照求得的网格尺度上小麦最接近于真实分布的概率进行分配，得到最接近真实分布情况的空间布局，即每个网格单元内的小麦种植比例或者种植面积。小麦先验分布参数pijl需先依据步骤S1得到的小麦种植适宜性面积数据以及由人口密度数据计算得到的市场可达性计算得到。

县域尺度小麦的播种面积、单产数据均来源于不同的农业统计年鉴。由于县区尺度范围太大，相当于搜集到每个县域上的平均单产和总播种面积，但是具体在更小的单元上播种面积怎么分配、单产是怎样的并不清楚。通过本发明步骤S3的模拟计算，可以将搜集到的县区尺度小麦单产、播种面积分配到网格尺度上，从而体现了每个县区内播种面积和单产的差异性。同时保证了每个县域内的网格单元上总产加和等于县域小麦总产，每个网格单元上的播种面积加和等于县区小麦播种面积总量，最后分配到网格单元上，对于优化种植分布更加精确。

步骤S4、获取每个网格单元内的小麦施肥强度数据。

获取网格尺度上小麦种植所需的氮肥、磷肥、钾肥等肥料的施用强度数据，即单位面积施肥量(单位：kg/ha)。

通过整理黄淮海区范围内的农业气象站点资料，获取各个农业气象站点上小麦的氮肥、磷肥、钾肥施用强度数据；其次将小麦种植分区与省界叠加，划分为各个亚区，统计各亚区范围内各农业气象站点上小麦施肥强度数据的平均值，并将其作为该亚区的小麦施肥强度数据；最后，通过将每个亚区内的各个网格单元赋值为该平均施肥强度，即完成空间化。施用强度数据为单位面积内的肥料施用量，只需要将该数值赋值到每个网格单元，不需要通过其他算法分配。

氮肥、磷肥、钾肥的施肥强度数据来源于各农业气象站点的统计资料。

步骤S5、获取网格单元内的耕地面积和灌溉面积。

网格单元内的耕地面积和灌溉面积，表示该网格单元内所有作物当前实际占用的耕地面积和实施灌溉部分的耕地面积。通过土地利用数据，计算得到网格单元内耕地所占的百分比，进而计算网格单元内的耕地面积；从全球灌溉地区数字地图v4得到种植区的各个网格单元内灌溉百分比数据，进而确定网格单元内灌溉所占的耕地面积。

步骤S6-1，对步骤S1-S5获得网格单元数据转化为点状数据。

在进行步骤S6的优化建模之前，还需要对前述步骤S1-S5获得的基础数据进行处理，处理过程如下。

利用ArcGIS的创建渔网工具，创建种植区范围固定距离的渔网点，如1km×1km。创建渔网工具为内嵌到ArcGIS软件中的应用工具，可以创建包含矩形网格所组成网络的要素类。创建渔网需要三组基本信息：渔网的空间范围、行数和列数以及旋转的角度。

在本发明中，利用此工具输入的空间范围为目标种植区(本发明为黄淮海区的范围，具体包括北京市、天津市、河北省、山东省、河南省、安徽省和江苏省)；行数和列数设置中，可以通过给定的网格长度和宽度确定，选择1km作为长和宽进行确定；旋转角度设置不考虑，选择不旋转。通过将三组基本信息输入，会输出覆盖研究区范围的一组渔网点，这些渔网点之间的距离均为1km。

渔网点表示将各个网格单元(矩形单元)转换为点的形式，方便优化建模处理。依据创建的渔网点数据提取各渔网点S1-S5得到的各类指标数据，即小麦在灌溉和雨养种植条件下的适宜种植面积、小麦生长发育所需的最小灌溉需水量、当前种植布局下的种植数据(灌溉和雨养小麦的种植面积和单产数据)、当前种植布局下氮肥/磷肥/钾肥的施用强度数据、耕地面积和灌溉面积，作为后续步骤S8优化过程所需的基础数据。

本步骤中，渔网点的作用主要是将面状的网格单元信息转换到点状，方便后续优化建模过程中数据的输入。

由于后面的优化建模步骤不能输入网格单元上的信息，但是可以输入每个点的信息。前述步骤S1-S5得到各指标数据为面状的网格数据，不能直接输入后续的优化建模过程，需要先将网格数据转换成点形式，即每个网格单元中心对应一个渔网点，将每个网格单元内的有效信息(S1-S5得到的各指标数据)提取到各个渔网点上。本步骤通过创建渔网工具提取的渔网点数据是覆盖研究区范围的一系列点状数据。通过将网格单元上的数据提取到点上，方便后续算法的进行。

上述步骤S6中提取的数据格式为.dbf格式(ArcGIS软件中数据存储时的文件格式)，为了便于后续建模处理，将步骤S6提取的.dbf格式数据转存为.xlsx格式，并利用Python，将数据进一步修正为方便GAMS软件读取的.gdx格式的文件。

步骤S6，以小麦的灌溉需水总量最小为优化目标，基于预定的约束条件和前述步骤得到的数据进行种植布局优化，得到优化后各个网格单元内小麦的种植面积。

本步骤S6中，通过构建目标函数和约束条件，结合上述步骤S1-S5提取的各类基础数据，利用GAMS软件，优化模拟，输出满足目标函数和约束条件下的小麦种植布局，即确定优化后各个网格单元内小麦的种植面积。

具体来说，在GAMS软件中新建工程文件，读取优化算法，选择Conopt4解决器，并依据自身数据大小，设定程序运行的内存限制，然后运行优化算法。优化算法包括目标函数和约束条件。

其中，目标函数表示了本专利的种植布局优化方法的优化目标，本发明中，布局优化方法的优化目标设定为保证优化后种植区小麦灌溉需水总量最小。这里，小麦灌溉需水总量为步骤S2中得到的每个网格单元内小麦灌溉需水量乘以对应的每个网格单元优化后的种植面积，得到每个网格单元内的灌溉需水总量，再将种植区范围内所有网格单元上的灌溉需水总量相加，即得到优化后种植区小麦灌溉需水总量。

本发明的实施例中，目标函数的实现是通过调用GAMS系统中“最小化”函数实现的，该函数的原理即通过比较满足不同约束条件下的各种情景，找到目标最小的情景。

在GAMS软件实现优化算法模型中，优化目标采用下述目标函数(6-1)实现：

IRD_j＝∑_io_ijl×CropArea_jl×Dirri_i×10,l＝irri (6-1)

其中，IRD_j表示目标种植区内作物j的灌溉需水总量；o_ijl表示优化后网格单元i内作物j在灌溉I或雨养R种植条件l下的种植比例；CropArea_jl表示作物j在目标种植区实际的灌溉I或雨养R种植条件l下的种植面积；Dirri_i表示网格单元i上的灌溉需水量；10为转换系数。

为达到上述优化目标，本发明的优化方法设置了多个约束条件，下面详细介绍本发明优化方法的约束条件，这些优化条件以约束函数的形式设置在优化模型中。

下述各个约束条件和目标函数之间的联系主要是通过种植比例o_ijl进行联系。优化过程中，通过各个约束条件，可以输出多种优化后的种植比例o_ijl，通过将各种o_ijl应用到目标函数(6-1)中，计算不同种植比例对应的灌溉需水总量。通过比较不同的灌溉需水总量，确定一种灌溉需水总量最小的方案，然后输出此情景下的o_ijl。约束函数如下所示：

∑_lo_ijl×CropArea_jl≤Lucc_i (6-2)

公式(6-2)约束了在每个网格单元i上，灌溉和雨养种植条件l下分配的小麦种植面积的加和小于等于该网格单元上的耕地面积(S5得到结果)。o_ijl表示优化后网格单元i上作物j在灌溉或雨养种植条件l下分配的种植比例；CropArea_jl表示作物j在目标种植区实际的灌溉或雨养种植条件l下的总种植面积；Lucc_i表示网格单元i内的耕地面积；

CropArea_jl×o_ijl≤Irri_i，l＝irri (6-3)

公式(6-3)约束了每个网格单元i上，灌溉条件下分配的小麦种植面积小于等于该网格单元内的灌溉面积(S5得到的灌溉面积数据)。Irri_l表示网格单元i内的灌溉面积；

CropArea_jl×o_ijl≤SuitArea_ijl (6-4)

公式(6-4)约束了每个网格单元i上，灌溉或雨养种植条件l下分配的小麦种植面积小于等于该网格内灌溉或雨养种植条件下的小麦适宜性种植面积(步骤S1获得的数据)。SuitArea_ijl表示网格单元i内作物j在灌溉或雨养种植条件l下的适宜性种植面积；

o_ijl≤S_ijl，l＝rainfed (6-5)

公式(6-5)约束了优化后每个网格单元i上分配的雨养小麦种植比例小于等于原布局上每个网格单元内的雨养小麦种植比例(步骤S3得到的当前布局下网格尺度的雨养小麦的种植比例)。s_ijl表示优化前网格单元i上作物j在种植条件l为雨养下的种植比例；rainfed表示种植条件l为雨养；

本发明的方法中，主要的优化目标是灌溉需水量最小。通过上述约束条件，使得优化后雨养小麦的种植比例不超过原布局下的雨养小麦种植比例。例如，黄淮海区江苏省和安徽省南部地区降雨充足，且以种植水稻为主。在实现灌溉水最小目标的背景下，如果不限制雨养小麦的种植比例，这些地区的雨养小麦会大幅增长，与实际的种植情况不符(以水稻种植为主，小麦种植面积较小)。因此设置约束，保证优化后雨养小麦的种植比例不增加。

公式(6-6)约束了优化后的小麦种植区保证在原小麦种植网格单元范围内(S3输出的当前小麦种植布局上的面积大于0的地区)。该约束保证了在网格单元i上，优化前作物j在灌溉I或雨养R种植条件l下分配的种植面积为0时，优化后该网格单元i上分配的同种种植条件l下的种植比例也为0；本约束条件限定了优化后小麦种植区域必须在原小麦种植布局范围内。当前小麦的种植布局是经过农户多年种植过程中所形成的，在当前小麦种植布局范围外，可能并不适合小麦的种植。通过约束保证了优化后小麦仍然在原有布局范围内，防止小麦在气候条件或人为条件限制下不适合的地方种植；

∑_i∈Z∑_lo_ijl×CropArea_jl×yld_ijl≥∑_i∈Z∑_lm×s_ijl×CropArea_jl×yld_ijl (6-7)

公式(6-7)约束了优化后小麦总产大于等于不同减产比例与优化前种植区小麦总产(S3步骤中，当前小麦单产和种植面积计算得到小麦总产量)的乘积。该约束中，分别设置了小麦总产不变、减产10％、20％和30％四个比例，分别用于约束优化后的总产量需求。一方面，产量不变条件即分析在保证现有小麦产量下的布局情况，另一方面，通过设置不同的减产百分比，分析在可接受的减产比例下的可能布局方案；

∑_i∈Z∑_lo_ijl×CropArea_jl×f_ij≤∑_i∈Z∑_ls_ijl×CropArea_jl×f_ij (6-8)

公式(6-8)约束了优化后种植区肥料(N、P、K)施用总量(步骤S4输出的肥料施用强度乘以优化后得到的小麦种植面积，得到优化后的肥料施用总量)小于等于优化前的总施用量(步骤S4输出的肥料施用强度结合S3得到的小麦种植面积，计算得到优化前肥料的总施用量)。该约束保证了优化后肥料的施用量不增加，减少对环境的影响；

∑_i∈Z∑_lo_ijl×CropArea_jl×Dirri_i≤∑_i∈Z∑_ls_ijl×CropArea_jl×Dirri_i (6-9)

公式(6-9)约束了优化后种植区小麦灌溉需水总量(步骤S2得到的小麦灌溉需水量乘以优化后的种植面积得到)小于等于优化前的小麦灌溉需水总量(S2得到的小麦灌溉需水量乘以优化前的小麦种植面积得到)，即优化后灌溉水总需求量不增加。通过该指标的约束，有利于达到本发明的优化目标，即保证优化后的灌溉需水量最小。通过限制优化后的水需求不增加，以实现灌溉水不超过原布局的最小目标。

步骤S7、优化后的结果进行空间化。

通过上述一系列约束算法进行优化，可以得到优化后各渔网点上的灌溉小麦种植面积和雨养小麦种植面积。步骤S6-1中通过创建的渔网点数据，提取每个渔网点对应的网格单元上的信息，即将面状的网格单元转换为了点状数据。本步骤S7中，通过ArcGIS软件的点转网格工具实现空间化，即通过点转栅格工具，将点数据重新转换为网格单元对应的面状数据形式，从而得到网格尺度下的优化结果。

本发明方法的最终目标是要实现网格尺度下的作物种植布局优化模拟，

但是由于优化建模过程中需要点状形式的数据输入，因此，本本发明先获取网格尺度的基础数据，然后转换为点状的，最后在转换为网格尺度的结果。

如上所述，本发明提出了一种网格化作物种植布局优化方法，以保证优化后种植区作物灌溉需水总量最小，用于缓解当前地区用水紧缺以及地下水位下降等生态环境问题。本发明的方法进一步考虑了不同种作物减产限制条件，分别模拟了作物总产不变、减少10％、减少20％以及减少30％情景下以节水为目标的可能布局。

本发明方法中的优化算法是基于网格单元进行的，通过非线性优化求解，以地区种植作物的灌溉需水总量最小为目标，同时设定了8个约束条件，分别保证(1)在每个网格单元上灌溉和雨养种植条件下分配的作物种植面积加和小于等于该网格单元上的耕地面积、(2)每个网格单元上分配的作物灌溉面积小于等于该网格单元上的灌溉面积、(3)每个网格单元上分配的作物灌溉面积或雨养面积小于等于网格单元上灌溉或雨养种植条件下的作物适宜性种植面积、(4)每个网格单元上分配的作物雨养面积应该小于等于原布局上的作物雨养面积、(5)优化后作物的种植区保证在原作物种植网格单元范围内、(6)优化后作物的总产大于等于不同缩减比例与优化前种植区小麦总产的乘积、(7)优化后种植区肥料总施用量小于等于优化前的总施用量、(8)优化后种植区作物灌溉需水总量小于等于优化前的作物灌溉需水总量。

通过上述约束条件的组合约束，能够保证得到更加合理的优化布局结果，并且能够在网格单元上精准详细的提出作物种植面积缩减目标，特别是精确到每个网格单元内的种植面积精准控制，从而达到节水目的。此外，该方法考虑了农民种植习惯和灌溉设施等实际情况，能够在不改变农户意愿的情景下调整作物种植布局，实施起来较为可行。

下面结合实施例对本发明的技术方案进行说明。

本发明中，以黄淮海区为例实施在不同产量目标下的网格化小麦种植布局优化方法。图3-5中分别显示了灌溉小麦、雨养小麦以及小麦总种植面积的空间分布，其中OPT表示优化后小麦种植面积的空间分布，CG表示优化后较优化前的种植面积变化量，INT表示优化前小麦种植面积的空间分布。0％、10％、20％和30％分别表示总产不变以及减产10％、20％和30％的情景。

1.种植区概况

黄淮海区(29°41′–42°37′N，110°21′–122°42′E)包含我国的第二大平原(黄淮海平原)，也是世界上最大的含水层之一。其面积大约31×10⁴km²，包含北京市、天津市、河北省、山东省、河南省、安徽省以及江苏省。从北到南，分别分布着海河、黄河和淮河流域。作为我国重要的农业区，黄淮海区对保障我国的粮食安全发挥着重要作用。该区域内，粮食作物总产量约占全国总产量的35％以上，而小麦的总产量甚至超过全国总产量的80％。黄淮海区主要作物种植制度为冬小麦-夏玉米，但由于黄淮海区受到季风气候的影响，年内降雨主要集中在7、8月份，冬小麦生育期内降雨匮乏且耗水量大，为保证小麦的稳定高产，当前多通过大量开采地下水灌溉来满足小麦的生长需求。然而，过度开采地下水引起了黄淮海区地下水位大幅下降和地下水漏斗群的产生，已有研究指出黄淮海区地下水位平均每年下降达0.8m，形成了总面积达4.1×10⁴km²的浅层地下水漏斗区，5.6×10⁴km²的深层地下水漏斗区，该现状严重制约地区内农业的可持续发展。基于黄淮海区面临的严重问题，我们选择黄淮海区作为实例进行解释本研究方法的可行性和有效性。

2.研究数据

(1)气象数据

研究使用1981-2015年期间黄淮海区平均温度、降雨、相对湿度、风速、大气压和日照时数等日值数据。其中，站点尺度气象数据来源于中国气象数据网，并使用ANUSPLIN软件对各项气象因子进行空间插值，插值结果空间分辨率为1km。网格尺度气象数据来源于国家青藏高原数据中心。

(2)作物统计观测数据

本研究使用的作物统计数据主要包括黄淮海区各农业气象站点的农作物生育期和产量观测数据；各类农业统计年鉴中获取1981年以来的黄淮海区县级尺度主要作物的单产、总产和播种面积数据。

(3)资源环境数据

研究使用的资源环境数据主要包含土壤数据、地形数据、土地利用数据、土地利用网格百分比数据。其中地形数据来源于地理空间数据云，空间分辨率为30m；土地利用数据来源于国家土地利用/土地覆被数据集(NLCD)；土壤类型、土壤pH、土壤有机碳等土壤属性数据来源于中国科学院南京土壤所的中国土壤数据集和中国科学院资源环境科学数据中心的土壤属性数据集，土壤有效含水量数据来源于ISRIC-WISE全球土壤属性数据集。

(4)管理和人口密度数据

本研究使用的管理数据包括黄淮海区化肥施用量数据、县级灌溉面积统计数据和灌溉面积网格数据。其中，作物施肥量数据和县级灌溉面积数据来源于农业统计年鉴，灌溉百分比数据来源于全球灌溉地区数据地图v4，分辨率为5′×5′；人口密度数据来源于世界人口分布图，空间分辨率为0.5′×0.5′。

(5)作物系数K_c值数据

本研究主要作物蒸散计算相关的作物系数K_c来源于《北方地区主要农作物灌溉用水定额》。该数据中包括黄淮海区范围主要作物的作物系数，该系数值以FAO的作物系数为标准，校正后与国内试验资料相对比，使得结果既可以验证FAO的标准方法在我国北方地区的实用性，也可使得结果具有更好的精度和更强的可比性。

(6)作物灌溉、雨养比例数据

本研究使用的灌溉、雨养比例数据来源于SPAM2010数据集。该数据集提供了2010年全球42种作物在灌溉条件和雨养种植条件下的种植面积空间分布，研究通过灌溉条件下作物种植面积与总种植面积的比例来确定黄淮海区的灌溉、雨养比例。

3.研究方法

(1)小麦潜在产量和适宜性面积模拟

结合各类气象数据、土壤数据、地形数据等，通过全球农业生态区划(GAEZ)模型，通过在气候条件限制、土壤水条件限制、病虫害和霜冻等灾害限制、土壤属性、地形、耕地比例和灌溉面积等共同限制下，模拟输出作物的生产潜力，同时输出的数据还包含灌溉和雨养种植条件下的适宜性种植面积数据。

(2)小麦最小灌溉需水量的计算

利用各类气象数据，采用FAO推荐的以能量平衡和水汽扩散理论为基础的Penman-Monteith公式，计算参考作物潜在蒸散量。并结合作物系数和水分胁迫系数，计算小麦在生育期内逐日的潜在蒸散量和实际蒸散量。其次结合降雨数据计算逐日有效降雨，小麦生育期内潜在蒸散量与有效降雨两者的差值即为小麦的最小灌溉需水量。具体算法参见步骤S2的描述。

(3)网格尺度上小麦种植面积和实际产量的模拟

小麦产量空间模拟方面，以最小化交叉信息熵和作物产量空间分布模型(SPAM)为基础，结合小麦的播种面积和单产等统计数据以及其他的气象、土壤、地形等数据模拟得到。作物产量空间分配模型首先需要将小麦的播种面积、单产等数据按照灌溉与否划分为灌溉和雨养种植条件下的对应数据；其次通过GAEZ模型模拟得到的各级小麦适宜性面积进行加和，得到最终适宜性面积数据；第三，通过人口密度等数据计算市场可达性；第四，基于作物的市场可达性、适宜性面积以及潜在单产等数据，计算网格单元尺度上的作物生产的成本参数和先验分布参数，越高的成本参数代表网格单元中分配的比例越大；第五，基于交叉信息熵理论，模拟得到网格单元上小麦最可能接近真实分布的概率；最后结合得到的接近真实的分布概率与种植面积数据、单产数据综合，计算得到网格尺度上的作物种植面积和单产。

(4)施肥数据空间化

通过统计年鉴和农业气象站点等统计数据，整理黄淮海区各省市及各农业气象站点的小麦氮肥、磷肥和钾肥的施用强度，借助ArcGIS工具，将统计数据进行空间化显示。

(5)小麦优化布局算法

为了得到小麦最小化灌溉用水需求的优化布局，分别模拟了小麦总产不变、减少10％、20％和30％的情况下，肥料用量不增加、网格尺度上雨养种植区面积不增加等条件下的优化布局。该方法使用了通用代数建模系统(GAMS)的Python API接口，选择Conopt4解决器中的非线性求解算法，模拟小麦的优化布局。

该方法首先通过构建渔网点，并提取各渔网点上的所需参数，如小麦的最小灌溉需水量、灌溉百分比网格数据、耕地百分比网格数据、优化前灌溉小麦和雨养小麦的分配比例、灌溉和雨养条件下的适宜性种植面积、灌溉和雨养条件下的小麦实际单产、灌溉和雨养条件下的小麦总种植面积、氮肥/磷肥/钾肥的施用量。最后，结合GAMS软件的非线性解决器进行优化求解，具体优化求解的目标函数及约束条件参见步骤S6。

(6)通过GAMS优化求解能够得到各渔网点上小麦种植面积的优化结果。通过ArcGIS的点转网格工具，可以将结果进行空间化，得到网格尺度上灌溉和雨养种植条件下的小麦优化布局结果。

4.研究结果

通过对黄淮海区小麦的种植布局进行优化模拟，输出满足最小化种植区灌溉用水需求的结果，保证小麦总产不变以及减少10％、20％和30％，同时施肥量不增加且小麦调减或调增均在小麦原网格布局上进行的优化布局，并对比分析了布局变化特征。

(1)黄淮海区优化布局的灌溉小麦种植面积空间分布变化

从种植区整体上看，在小麦总产不变以及分别缩减10％、20％和30％的情况下，灌溉小麦种植比例在灌溉需水量较大的中部地区大幅减少，大部分地区灌溉小麦种植比例减少30％以下，但是在山东省北部、河南省北部地区以及京津部分地区灌溉小麦缩减幅度在30％以上，同时河北省北部地区大部分灌溉小麦的种植面积缩减为0，但在部分小麦高产区灌溉比例有所增加；在降雨较为充沛、小麦灌溉需水量小的江苏省和安徽省南部，灌溉小麦比例有所增加。分省来看，京津地区灌溉小麦缩减幅度超过65％，调减比例在各省中最高。山东省、江苏省和安徽省的灌溉小麦缩减比例次之，分别可达各省市原布局比例的32.2-66％、29.5-30.9％和34.9-44.5％。河北省灌溉小麦缩减比例略低于山东省、江苏省和安徽省，在不同总产限制条件下降幅可达1122.13-1781.63×10³ha，在总产不变情景下甚至略有提升，增加250×10³ha。

(2)黄淮海区优化布局的雨养小麦种植面积空间分布变化

图4是根据本发明的布局优化方法实施后雨养小麦种植面积空间分布图。

如图4所示，黄淮海区雨养小麦的种植面积在不同减产百分比条件下的布局与原布局相差不大。其中江苏省雨养小麦减少面积占原小麦种植面积的比例最大，不同总产限制条件下下降比例为11.7-11.8％，其余各省份下降比例均低于6％。安徽省雨养小麦缩减幅度最大，但低于10.29×10³ha。山东省和河南省次之，调减幅度为4.96-7.43×10³ha，其余各省市调减幅度均低于1×10³ha。

(3)黄淮海区优化布局的小麦总种植面积空间分布变化

图5是根据本发明的布局优化方法实施后小麦总种植面积空间分布图。

如图5所示，在约束小麦总产不变以及分别减少10％、20％和30％的条件下，小麦的种植布局与原布局有较大差别，尤其是灌溉小麦的种植面积变化比较大。从不同省市的角度上看，京津地区总种植面积调减较少，不同总产限制条件下下降幅度均低于114×10³ha，但降幅占原京津地区小麦种植面积的60％以上，调减比例在各省市中最高。山东省、江苏省和安徽省小麦总种植面积调减幅度次之，分别可达各省市原布局的30.6-62.6％、26.3-27.4％、34.9-44.4％，其中山东省小麦总缩减面积最大，达到1211.93-2478.44×10³ha。河北省小麦总种植面积调减幅度略低于山东省、江苏省和安徽省，但随着总产不断下降，灌溉小麦的调减幅度增加最为明显。河南省虽然小麦种植面积调减比例略低于其他省市，但由于优化前小麦总面积就很高，因而调减幅度也处于较高水平，不同总产限制条件下可达1129.56-1788.3×10³ha；在总产不变约束下，河北省总播种面积可有提升。

(4)黄淮海区优化布局的灌溉需水总量变化

图6是根据本发明的布局优化方法实施后小麦灌溉需水量空间分布图。

如图6所示，从黄淮海区整体上看，不同总产限制条件下，大部分地区小麦灌溉需水量均有较大幅度下降，尤其是黄淮海区中部地区。其中，在河北省东部地区、山东省北部和河南省北部地区灌溉水消耗量减少超过10×10⁴m³。通过优化布局，在总产不变以及下降10％、20％和30％的条件下，黄淮海区整体可分别节约0.35-2.18×10¹⁰m³的灌溉用水。从不同省市来看，小麦产量较高但灌溉需水量低的江苏省由于灌溉小麦的面积有所增加，灌溉需水量在不同总产限制条件下有所增长，增加幅度为2.8×10⁸-1.27×10⁹m³。但由于该地区降雨相对充足，灌溉需水量少，增加的灌溉用水并不会增加黄淮海区整体的用水负担。山东省的小麦灌溉用水量缩减最大，不同总产限制条件下缩减幅度高达1.79-9.66×10⁹m³，缩减部分占山东省原灌溉需水量的13.4-72.2％。

5结论

通过限制小麦总产不变以及分别减少10％、20％和30％，对黄淮海区小麦的种植布局进行优化调整。在不同减产百分比条件下，小麦种植比例在灌溉需水量较大的中部地区调减明显，全黄淮海区可分别缩减0.35-2.18×10¹⁰m³的灌溉用水。在总产减少30％的情景下节水量最大，但是节约单位体积的灌溉用水所需付出的总产减产代价最高；而在总产限制减少10％的条件下，虽然节水量少，但是每节约1m³的灌溉用水需付出的减产代价为1.11kg/m³。经过优化布局调整后的大部分小麦种植区内平均水分利用效率提升，且区域内产量差有所下降。因而，对高耗水的山东省、河南省和河北省等地，建议对小麦进行适当的流出，从而缓解区域用水压力，促进黄淮海区水资源的高效利用。

本领域普通技术人员可以理解实现上述实施例方法中的全部或部分步骤是可以通过程序来指令相关硬件来完成，所述的程序可以存储于一计算机可读取存储介质中，该程序在执行时，包括如上述各方法的实施例的流程。其中，所述的存储介质可为磁碟、光盘、只读存储记忆体(Read-OnlyMemory，简称ROM)或随机存取存储器(RandomAccessMemory，简称RAM)等。

应当理解的是，本发明的上述具体实施方式仅仅用于示例性说明或解释本发明的原理，而不构成对本发明的限制。因此，在不偏离本发明的精神和范围的情况下所做的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。此外，本发明所附权利要求旨在涵盖落入所附权利要求范围和边界、或者这种范围和边界的等同形式内的全部变化和修改例。

Claims (13)

1.一种网格化作物种植布局优化方法，包括：

步骤S1、对目标种植区域建立网格模型，获取每个网格单元内的作物潜在单产和适宜种植面积；

步骤S2、获取每个网格单元内的作物最小灌溉需水量；

步骤S3、根据目标种植区的县级统计种植数据，分配得到每个网格单元内的种植数据，种植数据包括灌溉和雨养种植条件下的作物种植面积和单产数据；

步骤S4、获取每个网格单元内的作物施肥强度数据；

步骤S5、获取网格单元内的耕地面积和灌溉面积；

步骤S6，以灌溉需水总量最小为优化目标，基于预定的约束条件和前述步骤S1-S5得到的数据进行种植布局优化，得到优化后各个网格单元内的作物种植面积；

其中，每个网格单元内作物最小灌溉需水量乘以对应的优化后每个网格单元内的种植面积，得到每个网格单元内的灌溉需水总量，再将种植区范围内所有网格单元上的灌溉需水总量相加，即得到优化后种植区作物灌溉需水总量；

其中，通过下述约束条件实现所述优化目标：

(1)在每个网格单元上灌溉和雨养种植条件下分配的作物种植面积加和小于等于该网格单元上的耕地面积；

(2)每个网格单元上分配的作物灌溉面积小于等于该网格单元上的灌溉面积；

(3)每个网格单元上分配的作物灌溉或雨养面积小于等于网格单元上灌溉或雨养种植条件下的该作物适宜性种植面积；

(4)每个网格单元上分配的作物雨养面积小于等于原布局上该网格单元内的作物雨养面积；

(5)优化后作物的种植区在原作物种植网格单元范围内；

(6)优化后作物的总产大于等于不同减产比例与优化前种植区作物总产的乘积；

(7)优化后种植区肥料总施用量小于等于优化前的总施用量；

(8)优化后种植区作物灌溉需水总量小于等于优化前的作物灌溉需水总量。

2.根据权利要求1所述的方法，所述步骤S2进一步包括：

步骤S21，计算网格单元内的参考作物潜在蒸散量；

步骤S22，结合参考作物潜在蒸散量、作物系数和水分胁迫系数，计算生育期内作物潜在蒸散量和作物实际蒸散量，其中作物潜在蒸散量基于参考作物潜在蒸散量和作物系数计算得到，作物实际蒸散量则基于作物潜在蒸散量和水分胁迫系数计算得到；

步骤S23，计算网格单元内的有效降雨量；

步骤S24，基于作物潜在蒸散量和有效降雨量计算每个网格单元内的作物最小灌溉需水量。

3.根据权利要求1所述的方法，所述步骤S3进一步包括：

步骤S31，计算灌溉和雨养种植条件下作物的单产和播种面积；

步骤S32，计算作物在灌溉和雨养种植条件下的适宜种植面积；

步骤S33，计算作物的市场可达性；

步骤S34，计算作物先验分布参数；

步骤S35，基于最小化交叉信息熵方法求解作物在网格单元上分配的种植比例；

步骤S36，计算每个网格单元上作物的单产和种植面积。

4.根据权利要求1所述的方法，在所述步骤S6之前，还包括步骤S6-1，对步骤S1-S5获得的网格数据转化为点状数据。

5.根据权利要求1所述的方法，所述优化目标进一步表征为下述目标函数：

IRD_j＝∑_io_ijl×CropArea_jl×Dirri_i×10，l＝irri (6-1)

其中，IRD_j表示目标种植区内作物j的灌溉需水总量；o_ijl表示优化后网格单元i内作物j在灌溉或雨养种植条件l下的种植比例；CropArea_jl表示作物j在目标种植区实际的灌溉或雨养种植条件l下的总种植面积；Dirri_i表示网格单元i上的最小灌溉需水量；10为转换系数；irri表示种植条件，l为灌溉和雨养种植条件。

6.根据权利要求1所述的方法，所述8个约束条件分别采用下述约束函数实现：

∑_lo_ijl×CropArea_jl≤Lucc_i (8-1)

公式(8-1)约束了在每个网格单元i上，灌溉和雨养种植条件l下分配的作物j种植面积加和小于等于该网格单元i上的耕地面积，o_ijl表示优化后网格单元i上作物j在灌溉I或雨养R种植条件l下分配的种植比例，CropArea_jl表示作物j在目标种植区实际的灌溉I或雨养R种植条件l下的总种植面积，Lucc_i表示网格单元i内的耕地面积；

CropArea_jl×o_ijl≤Irri_i，l＝irri (8-2)

公式(8-2)约束了每个网格单元i上，灌溉种植条件下l＝irri分配的作物j种植面积小于等于该网格单元i内的灌溉面积；Irri_i表示网格单元i内的灌溉面积；

CropArea_jl×o_ijl≤SuitArea_ijl (8-3)

公式(8-3)约束了每个网格单元i上，灌溉或雨养种植条件l下分配的作物j种植面积小于等于该网格单元i内灌溉或雨养种植条件l下的作物j适宜性种植面积；SuitArea_ijl表示网格单元i内作物j在灌溉或雨养种植条件l下的适宜性种植面积；

o_ijl≤s_ijl，l＝rainfed (8-4)

公式(8-4)约束了优化后每个网格单元i上分配的作物j雨养面积小于等于原布局上每个网格单元i内的作物j雨养面积；s_ijl表示优化前网格单元i上作物j在种植条件l为雨养情况下的种植比例；rainfed表示种植条件l为雨养；

公式(8-5)约束了优化后的作物种植区保证在原作物种植网格单元范围内，即在网格单元i上，优化前作物j在灌溉或雨养种植条件l下分配的种植面积为0时，优化后该网格单元上分配的同种植条件下的种植比例也为0；

∑_i∈Z∑_lo_ijl×CropArea_jl×yld_ijl≥∑_i∈Z∑_lm×s_ijl×CropArea_jl×yld_ijl (8-6)

公式(8-6)约束了优化后作物j总产大于等于不同减产比例与优化前种植区作物j总产的乘积；yld_ijl表示网格单元i上作物j在灌溉或雨养种植条件l下的单产；m表示减产比例；Z表示研究区范围；

∑_i∈Z∑_lo_ijl×CropArea_jl×f_ij≤∑_i∈Z∑_ls_ijl×CropArea_jl×f_ij (8-7)

公式(8-7)约束了优化后种植区作物j的肥料施用总量小于等于优化前作物j的肥料施用总量；f_ij表示网格单元i内作物j的施肥强度数据；

∑_i∈Z∑_lo_ijl×CropArea_jl×Dirri_i≤∑_i∈Z∑_ls_ijl×CropArea_jl×Dirri_i (8-8)

公式(8-8)约束了优化后种植区作物j的灌溉需水总量不超过优化前的作物j灌溉需水总量，Dirrii表示网格单元i内的最小灌溉需水量。

7.根据权利要求2所述的方法，所述步骤S21中，采用下述公式计算网格单元内的参考作物潜在蒸散量ET₀：

式中，Δ表示为温度曲线和饱和水气压相关的斜率，KPa/℃；R_n表示地表净辐射量，MJ/m²·d；G表示土壤热通量，MJ/m²·d，一般取0；γ表示干湿度常数，KPa/℃；f表示空气平均温度，℃；U₂表示2m处风速，m/s；e_s表示饱和水气压，KPa；e_a表示实际观测水气压差，KPa。

8.根据权利要求7所述的方法，所述步骤S22中，采用下述公式计算作物潜在蒸散量PET_c，t和作物实际蒸散量ET_a，t：

PET_c，t＝K_c，tET_0，t (10-1)

ET_a，t＝K_c，tK_s，tET_0，t (10-2)

公式(10-1)为作物在第t天潜在蒸散量PET_c，t的计算，其中K_c，t为第t天的作物系数，ET_0，t为第t天的参考作物潜在蒸散；

公式(10-2)为作物在第t天实际蒸散量ET_a，t的计算，其中K_s，t为第t天的作物水分胁迫系数，具体计算公式如下：

式中，S_t为作物根部平均深度的土壤含水量；S_max为根区最大土壤有效含水量；P为作物从根区土壤中吸收的无水分胁迫的水分含量占根区最大土壤有效含水量S_max的比例；

S_t＝S_t-1+Δt(P_eff，t+I_t-ET_a，t-R_t) (10-4)

式中，S_t-1为前一时期的土壤含水量，Δt为1天，P_eff，t为第t天的有效降雨，I_t为第t天的灌溉用水；ET_a，t为第t天的作物实际蒸散量；R_t为第t天的径流；

R_t＝(P_eff，t+I_t)×(S_t/S_max)^γ (10-5)

式中，γ为常数，在雨养地区为2，灌溉地区为3；

P＝P_std+0.04(5-PEr_c，t) (10-6)

式中，P_std为常数，PET_c，t为作物在第t天的潜在蒸散量。

9.根据权利要求8所述的方法，所述步骤S23中，采用下述公式计算网格单元内的有效降雨量：

P_eff，t＝P_t-ΔSnow (11-1)

式中，P_t表示第t天的实际降雨量，当第t天的平均温度低于0℃时，认为ΔSnow等于P_t，当第t-1天有积雪，且第t天的平均温度高于0℃，则认为ΔSnow等于平均温度乘以4mm/℃。

10.根据权利要求9所述的方法，所述步骤S24中，采用下述公式计算网格区域的作物最小灌溉需水量Dirri：

当生育期内的作物潜在蒸散量∑PET_c，t大于等于有效降雨量∑P_eff，t时，最小灌溉需水量Dirri等于两者的差值，当生育期内作物潜在蒸散量∑PET_c，t小于有效降雨量∑P_eff，t时，则灌溉需水量Dirri为0。

11.根据权利要求3所述的方法，所述步骤S31中，通过下述公式计算县级尺度灌溉和雨养种植条件下作物的单产和播种面积：

Yield_jI＝∝_j×Yield_jR (13-2)

其中，Yield_j为作物j的县级统计单产；Yield_jR和Yield_jI分别表示各县域作物j在雨养R和灌溉I种植条件l下的单产水平；Percent_jR和Percent_jI分别表示各县域作物j在雨养R和灌溉I种植条件l下的种植比例，即雨养作物种植面积和灌溉作物种植面积占总种植面积的比例；∝_j为用于计算作物j灌溉单产和雨养单产的比例参数；

HarvestArea_jl＝Percent_jl×HA_j (13-3)

其中，HarvestArea_jl表示县级尺度作物j在灌溉或雨养种植条件l下的播种面积，HA_j表示县级统计数据中作物j的总播种面积；Percent_jl表示作物j在灌溉或雨养种植条件l下的种植比例。

12.根据权利要求3所述的方法，所述步骤S35中，通过下述公式基于最小化交叉信息熵方法求解网格单元上的作物分配比例：

En_ijl＝∑_i∑_j∑_ls_ijllns_ijl-∑_i∑_j∑_ls_ijllnp_ijl (14-1)

其中，En_ijl为交叉熵，s_ijl表示网格单元i上分配的作物j在不同种植条件l下的种植比例；p_ijl表示网格单元i上作物j在不同种植条件l下分配的初始种植比例，即作物先验分布参数；

在进行最小交叉熵的模拟时，需要同时满足以下约束条件：

∑_is_ijl＝1 (14-2)

∑_lCropArea_jl×s_ijl≤lucc_i (14-3)

CropArea_jl×s_ijl≤SuitArea_ijl (14-4)

∑_i∈c∑_lCropArea_jl×s_ijl＝StArea_jc (14-5)

∑_lCropArea_jl×s_ijl≤IrriArea_i，l＝irri (14-6)

其中，IrriArea_i为网格单元i上的灌溉面积；lucc_i为网格单元i上的耕地面积；CropArea_jl表示目标种植区内作物j在灌溉或雨养种植条件l下的总种植面积；SuitArea_ijl表示网格单元i上作物j在灌溉或雨养种植条件l下的适宜性种植面积；StArea_jc表示县域c统计数据中作物j实际播种面积与作物j熟制的比例；公式(14-2)约束了作物j在灌溉和雨养种植条件l下，各县域内所有网格单元i上分配比例的加和等于1；公式(14-3)约束了作物j在网格单元i上，在灌溉和雨养种植条件l下分配的作物种植面积加和小于等于网格单元i内的耕地面积；公式(14-4)约束了作物j在网格单元i上不同灌溉或雨养种植条件l下分配的作物种植面积小于等于作物适宜种植面积；公式(14-5)约束了各县域内作物j在所有网格单元i上分配的种植面积等于统计数据中作物j实际播种面积与作物j熟制的比例；公式(14-6)约束作物j在网格单元i上灌溉种植条件下分配的作物种植面积小于等于该网格单元i上的灌溉面积；

其中，CropArea_jl的计算公式如下：

其中，CropIntensity_jl为作物j在灌溉或雨养种植条件l下的熟制。

13.根据权利要求12所述的方法，所述步骤S36中，通过下述公式计算每个网格单元上作物的单产和种植面积：

依据求得的s_ijl，可以求得各网格单元上分配的种植面积Area_ijl；

Area_ijl＝s_ijl×CropArea_jl (15-1)

依据作物j在灌溉或雨养种植条件l下的潜在单产PotYield_ijl，可以求得各县域内的平均潜在单产Yavg_jl，并依据公式(15-3)计算得到各网格单元i上作物j在不同灌溉或雨养种植条件l下的实际单产Y_ijl；

其中，Yield_jl表示县级尺度上作物j在灌溉或雨养种植条件l下的单产水平。

Images