一种土地布局优化方法

1.本发明属于土地布局优化技术领域，具体涉及一种土地布局优化方法。


背景技术：

2.随着全球粮食需求的增加，耕地持续扩张，占用原本森林、草原及湿地等生态系统，破坏了野生动植物的栖息地和迁徙路径，成为全球自然栖息地和生物多样性丧失的主要驱动因素。而随着全球城市人口爆炸式增长，对城市用地的需求也在不断增加，并成为当前全球土地利用布局变化的主要驱动因素，也是直接导致生物栖息地退化和破碎、生物多样性丧失的主要原因之一。此外，不断扩张的耕地和城市用地还会进一步的造成生物量和土壤碳损失，导致碳储量服务减少。
3.在已有的以实现未来粮食安全或生态安全的土地利用布局模拟研究中，往往仅通过对耕地这一单一类型进行分配以实现粮食安全或生态安全，而忽略了建设用地与耕地或生态用地等多类型土地利用之间的权衡协同关系。其次，在多类型的土地利用布局模拟研究中，通常是在模拟未来土地利用布局基础上，评估多类型土地利用布局对粮食生产或者生态服务的影响程度，并不能保证未来的土地利用布局能够实现粮食安全和生态安全。
4.因此，如何平衡耕地、建设用地以及生态用地等多类型土地利用之间的冲突，模拟未来最优的土地利用布局，在保障未来粮食安全和建设用地需求的基础上，最大程度保证生物多样性或碳储量等生态安全成为当前全球关注的重点问题。


技术实现要素：

5.本发明的目的是提供一种土地布局优化方法，通过设置目标函数和约束条件对目标土地布局进行优化，在实现作物总热量最大、总物种丰富度最大和总碳储量最大的同时，保证粮食安全和生态安全。
6.为解决上述问题，本发明提供了一种土地布局优化方法，包括：
7.设置优化方法的目标函数，所述目标函数的优化目标至少包括：实现作物总热量最大、总物种丰富度最大和总碳储量最大；
8.设置优化方法的约束条件，所述约束条件至少包括：优化后目标区域内每种类型作物的总热量不减少、优化后目标区域内总物种丰富度不下降，以及优化后目标区域内总碳储量不减少；
9.基于所述约束条件和目标函数，对目标土地布局进行优化。
10.优选地，约束条件还包括：
11.优化后目标区域内每种类型作物的种植面积不超过该类型作物的种植适宜性面积；
12.优化后目标区域内的耕地面积不低于目标耕地面积；
13.优化后目标区域内的建设用地面积占比不低于目标建设用地面积占比；
14.优化后目标区域内的自然保护区的总物种丰富度不减少；
15.优化后目标区域内的水域面积不变；
16.优化前未利用地面积占比超过80％或者海拔超过4500m的目标区域，优化后该目标区域内未利用地的面积不变。
17.优选地，目标函数的输入参数包括：
18.目标区域内各种土地利用类型的用地面积；
19.目标区域内各种类型作物的作物播种面积和作物单产；
20.目标区域内作物种植面积占用耕地的比例；
21.目标区域内各种土地利用类型对应的单位面积碳储量；
22.目标区域内土地利用类型转变后对应的单位面积物种丰富度；
23.目标区域内目标耕地面积和目标建设用地面积占比。
24.优选地，目标区域内各种土地利用类型的用地面积，通过下述步骤获取：
25.获取目标区域的土地利用百分比栅格数据；
26.利用所述土地利用百分比栅格数据，计算目标区域上每个栅格单元内各种土地利用类型的用地面积。
27.优选地，目标区域内各种类型作物的作物播种面积和作物单产，通过下述步骤获取：
28.获取目标区域的各种类型作物的作物产量数据集；
29.计算目标区域上某年份每个栅格单元内各种类型作物的作物播种面积；
30.计算目标区域上某年份每个栅格单元内各种类型作物的作物单产；
31.作物产量数据集包括：第一作物产量数据集和第二作物产量数据集；
32.第一作物产量数据集为具有各种作物类型的作物播种面积和作物单产的统计数据集；
33.第二作物产量数据集为具有各种作物类型的栅格尺度的作物播种面积和作物单产的数据集。
34.优选地，作物播种面积的分配算法为：
[0035][0036]
其中，r为栅格单元的编号；c为目标区域；
[0037]
f为作物类型；y为模拟年份；
[0038]
ha
r，c，f，y
表示模拟年份y中，目标区域c上的栅格单元r上作物f的播种面积；
[0039]
ha
r，c，f，2010
表示第二作物产量数据集中，目标区域c上栅格单元r上作物f的播种面积；
[0040]
staha
c，f，y
表示第一作物产量数据集中，模拟年份y中，目标区域c上栅格单元r上作物f的播种面积。
[0041]
优选地，作物单产的分配算法为：
[0042][0043]
[0044]
其中，y
r，c，f，2010
表示第二作物产量数据集中目标区域c上栅格单元r上作物f的单产；
[0045]
yavg
c，f，2010
表示由第二作物产量数据集计算得到的目标区域c上作物f的平均单产；
[0046]
stay
c，f，y
表示第一作物产量数据集中模拟年份y中，目标区域c上作物f的平均单产；
[0047]yr，c，f，y
表示模拟年份y中，目标区域c的栅格单元r上作物f的单产。
[0048]
优选地，各种土地利用类型对应的单位面积碳储量，通过下述步骤获取：
[0049]
在所述目标区域内设置采样点，获取各个采样点处的各种土地利用类型的单位面积碳储量数据；
[0050]
根据各采样点处的各种土地利用类型的单位面积碳储量数据，利用统计学方法获得栅格尺度各种土地利用类型的单位面积碳储量数据。
[0051]
优选地，目标区域内土地利用类型转变后对应的单位面积物种丰富度，通过下述步骤获取：
[0052]
获取基于当前年份土地分布情况下栅格尺度的单位面积总物种丰富度数据，并作为单位面积物种丰富度初始数值；
[0053]
计算各种土地利用类型之间转换带来的单位面积物种丰富度变化百分比；
[0054]
基于所述单位面积物种丰富度初始数值和所述单位面积物种丰富度变化百分比，计算得到变化后的目标区域内土地利用类型转变后对应的单位面积物种丰富度。
[0055]
优选地，目标区域的目标耕地面积和目标建设用地面积，通过下述步骤获取：
[0056]
获取目标区域内当前耕地面积占全国当前总耕地面积的比例；
[0057]
根据上述比例和全国目标总耕地面积获取目标区域的目标耕地面积；
[0058]
获取目标区域在不同时段的建设用地占土地利用总面积的比例；
[0059]
将不同时段的建设用地占比的平均值作为目标区域的目标建设用地面积占比。
[0060]
优选地，目标函数表征为：
[0061][0062]
其中，i表示不同的栅格单元，n表示目标区域栅格单元的总数，i＝1
…
n；
[0063]
j表示不同土地利用类型，j＝1
…
6；
[0064]
l表示作物类型；
[0065]
w1表示作物总热量的权重；
[0066]
w2表示总碳储量的权重；
[0067]
w3表示总物种丰富度的权重，w1+w2+w3＝1；
[0068]yi，l
表示栅格单元i上作物l的单产；
[0069]
cal
i，l
表示栅格单元i上每100g作物l可食用部分所含卡路里；
[0070]mi，l
表示栅格单元i上作物l的种植面积占耕地面积的比例；
[0071]oi，j
表示优化后栅格单元i上第j类土地类型的用地面积；
[0072]
x
i，j
表示优化前栅格单元i上第j类土地类型的用地面积；
[0073]ci，j
表示栅格单元i上第j类土地类型对应的单位面积碳储量；
[0074]bi，j
表示栅格单元i上第j类土地类型对应的单位面积物种丰富度。
[0075]
本发明的上述技术方案具有如下有益的技术效果：
[0076]
本发明的土地布局优化方法在栅格尺度上通过优化不同土地利用类型分布，通过约束作物的总热量不减少、物种丰富度不下降以及碳储量不减少以保证粮食安全目标和生态安全目标，同时约束未来建设用地分布、未来建设用地需求和未来耕地需求等，以实现对耕地、建设用地和生态用地等多地类进行综合调整。
附图说明
[0077]
图1是本发明的土地布局优化方法的流程图；
[0078]
图2是本发明目标函数设置方法的流程图；
[0079]
图3是利用软件优化建模流程的示意图。
具体实施方式
[0080]
为使本发明的目的、技术方案和优点更加清楚明了，下面结合具体实施方式并参照附图，对本发明进一步详细说明。应该理解，这些描述只是示例性的，而并非要限制本发明的范围。此外，在以下说明中，省略了对公知结构和技术的描述，以避免不必要地混淆本发明的概念。
[0081]
在本发明的描述中，需要说明的是，术语“第一”、“第二”、“第三”仅用于描述目的，而不能理解为指示或暗示相对重要性。
[0082]
当前的土地布局可以采用多种优化方法，例如：
[0083]
针对实现多目标的土地利用布局，首先确定不同区域的生态安全目标，然后在该区域内将某一种土地利用类型(一般为耕地)在适宜分布的地区进行优化分配，在实现生态安全目标的同时实现粮食安全目标。该技术体现为对多目标、单类型的优化，即虽然针对的目标是实现粮食安全和生态安全，但只体现在对某单一土地利用类型(如耕地)进行优化。
[0084]
另外，针对多类型土地利用的布局，通常做法是对多土地利用类型的预测。即，通过历史阶段两期或以上的土地利用布局数据，应用机器学习等算法进行综合比较，确定不同土地利用类型转化的概率，进而预测未来多类型土地利用布局方案。
[0085]
可见，这种针对多类型土地利用并对应多目标的研究，主要体现为对未来的影响。即，通过土地利用预测模型模拟未来多类型土地利用布局，进而评估未来土地利用布局下对生态系统服务或者粮食安全的影响。
[0086]
然而，上述方式具有一定的缺陷。
[0087]
例如，在以粮食安全或生态安全为目标的土地利用布局模拟中，主要考虑对耕地这一单一地类进行优化分配，以保证粮食安全和生态安全目标同时实现，而忽略了这种优化布局下，是否能够满足城市人口增加带来的建设用地需求。
[0088]
又例如，多类型土地利用布局模拟中，通过应用算法模拟历史时期不同地区土地利用布局转换的概率，进行未来多类型土地利用布局预测。这种方法虽然能够实现多类型的土地利用布局模拟，但是该算法是以土地利用模拟为基础的，只能分析未来多类型土地利用变化对粮食生产和生态服务等的影响程度，而不能实现未来的粮食安全和生态安全目标。
[0089]
可见，针对多类型土地利用并对应多目标的研究，只能评估影响，并不能实现同时达到粮食安全和生态安全多目标的土地利用类型的优化。
[0090]
基于上述问题和缺陷，本发明提出了一种土地布局优化方法，在栅格尺度上通过优化不同土地利用类型分布，同时实现未来粮食安全和生态安全目标。该方法通过约束作物的总热量不减少、物种丰富度不下降以及碳储量不减少以保证粮食安全目标和生态安全目标，同时约束未来建设用地分布、未来建设用地需求和未来耕地需求等，以实现对耕地、建设用地和生态用地等多地类进行综合调整。该算法一方面通过约束条件的限制，保证了三者的基本安全，又通过目标函数的设定，保证三者同时达到最优水平。有助于国家或地方政府更加方便的展开有针对性的、具体的国土空间规划措施，从栅格单元进行管理、调整，从而有效的协调耕地、建设用地和生态用地之间的竞争，以实现粮食安全和生态安全。
[0091]
实施例一：
[0092]
图1示出了本发明的一种土地布局优化方法的流程。
[0093]
如图1所示，该方法包括以下步骤：
[0094]
步骤s1、设置优化方法的目标函数，目标函数的优化目标包括：实现作物总热量最大、总物种丰富度最大和总碳储量最大；
[0095]
步骤s2、设置优化方法的约束条件，约束条件包括：优化后目标区域内每种类型作物的总热量不减少、优化后目标区域内总物种丰富度不下降，以及优化后目标区域内总碳储量不减少；
[0096]
步骤s3、基于目标函数和约束条件，对目标土地布局进行优化。
[0097]
本发明的土地布局优化方法，针对目标区域，设置约束条件以及目标函数对土地布局进行优化。在设定的目标函数和约束条件下，将收集的各类数据降尺度到栅格尺度，对每个栅格单元内的不同土地利用类型面积进行优化，得到优化后的栅格尺度的目标土地利用类型面积数据。在实现对目标土地布局优化调整的同时，保证粮食安全和生态安全。
[0098]
为了更加清楚地说明，下面对上述步骤进行详细地介绍。
[0099]
步骤s1、设置优化方法的目标函数；
[0100]
一个优选的实施例中，目标函数的输入参数包括：
[0101]
目标区域内各种土地利用类型的用地面积；
[0102]
目标区域内各种类型作物的作物播种面积和作物单产；
[0103]
目标区域内作物种植面积占用耕地的比例；
[0104]
目标区域内各种土地利用类型对应的单位面积碳储量；
[0105]
目标区域内土地利用类型转变后对应的单位面积物种丰富度。
[0106]
进一步的，步骤s1还包括获取上述输入参数的步骤，包括：
[0107]
步骤s11、计算目标区域内各种土地利用类型的用地面积；
[0108]
步骤s12、计算目标区域内各种类型作物的作物播种面积和作物单产；
[0109]
步骤s13、计算目标区域内作物种植面积占用耕地的比例；
[0110]
步骤s14、计算目标区域内各种土地利用类型对应的单位面积碳储量；
[0111]
步骤s15、计算目标区域内土地利用类型转变后对应的单位面积物种丰富度。
[0112]
一个优选实施方式中，步骤s11、目标区域内各种土地利用类型的用地面积，包括：
[0113]
步骤s111、获取目标区域的土地利用百分比栅格数据；
[0114]
步骤s112、利用土地利用百分比栅格数据，计算目标区域上每个栅格单元内各种土地利用类型的用地面积；
[0115]
其中，土地利用类型包括：耕地、林地、草地、建设用地、水域和未利用地，每个栅格单元内均可包含上述六类土地利用类型。
[0116]
土地利用百分比栅格数据，为地理信息系统中常用的数据结构，表示每个栅格单元内不同土地利用类型的面积占该栅格单元总面积的比例。
[0117]
栅格单元，即规则的网格单元，具有一定的尺寸和数值，尺寸对应一定的地理空间范围；栅格单元的大小表示栅格数据的空间分辨率。
[0118]
具体的，将目标区域划分成若干个空间分辨率为10km的栅格单元，每个栅格单元总面积为100km2。再结合土地利用百分比栅格数据，可计算得到每个栅格单元内各种土地利用类型的占地面积。将其作为第一项输入参数，用于保证粮食安全。
[0119]
一个优选实施方式中，步骤s12、计算目标区域内各种类型作物的作物播种面积和作物单产，包括：
[0120]
步骤s121、基于目标区域内各种类型作物的作物产量数据集，计算各目标区域上某年度每个栅格单元内各种不同类型作物的作物播种面积；
[0121]
步骤s122、计算目标区域上每个栅格单元内各种不同类型作物的作物单产；其中，作物产量数据集包括：第一作物产量数据集和第二作物产量数据集。
[0122]
具体的，第一作物产量数据集包括：各种作物类型的作物播种面积和作物单产数据，可以为源于各类农业统计年鉴的不同作物类型的作物播种面积和作物单产的统计数据。
[0123]
第二作物产量数据集包括：具有不同作物类型的栅格尺度的作物播种面积和作物单产数据，可以采用国际上公开发布的spam2010数据集。该数据集的分辨率约为10km
×
10km，是将数据从较粗的数据单位分配到栅格单元等较细的数据单位，揭示了作物表现的空间模式。
[0124]
其中，spam2010数据集包括：粮食作物、经济作物、水果等共41种作物。
[0125]
进一步的，步骤s121中，作物播种面积的分配算法为：
[0126][0127]
其中，r为栅格单元的编号；c为目标区域；
[0128]
f为作物类型；y为模拟年份；
[0129]
ha
r，c，f，y
表示模拟年份y中，目标区域c上的栅格单元r上作物f的播种面积；
[0130]
ha
r，c，f，2010
表示第二作物产量数据集中，目标区域c上栅格单元r上作物f的播种面积；
[0131]
staha
c，f，y
表示第一作物产量数据集中，模拟年份y中，目标区域c上栅格单元r上作物f的播种面积。
[0132]
步骤s122中，作物单产的分配算法为：
[0133]
[0134][0135]
其中，y
r，c，f，2010
表示第二作物产量数据集中目标区域c上栅格单元r上作物f的单产；
[0136]
yavg
c，f，2010
表示由第二作物产量数据集计算得到的目标区域c上作物f的平均单产；
[0137]
stay
c，f，y
表示第一作物产量数据集中模拟年份y中，目标区域c上作物f的平均单产；
[0138]yr，c，f，y
表示模拟年份y中，目标区域c的栅格单元r上作物f的单产。
[0139]
本步骤中，通过结合spam2010数据集中的栅格尺度的播种面积数据和单产数据，实现对该年度由大于栅格尺度的作物数据，到栅格尺度数据的降尺度，使最后的土地布局优化更加精确。本步骤中将播种面积和作物单产作为第二项输入参数，用于保证粮食安全。
[0140]
步骤s13、计算目标区域上作物种植面积占用耕地的比例，包括：
[0141]
步骤s131、计算目标区域上每个栅格单元内作物的种植面积；
[0142]
步骤s132、计算每个栅格单元内作物的种植面积与耕地用地面积的比例。
[0143]
具体的，步骤s131包括：
[0144]
步骤s1311、结合步骤s12中得到的目标区域每个栅格单元内各种类型作物的作物播种面积，计算每个栅格单元内栅格尺度的所有类型作物的播种面积；
[0145]
步骤s1312、基于搜集现有研究数据集中栅格尺度的作物熟制数据，计算每个栅格单元内栅格尺度的所有作物的种植面积；
[0146]
具体的，播种面积为考虑作物熟制后的作物总面积，通过将播种面积与作物熟制数据相除，即可得到作物的种植面积。
[0147]
步骤s132、由步骤s131得到每个栅格单元内栅格尺度的作物的种植面积之后，在通过与步骤s11中计算得到的每个栅格单元内土地利用类型为耕地的面积相除，即可得到每个栅格单元内作物的种植面积与耕地面积的比例，作为第三项输入参数，用于确保粮食安全。
[0148]
步骤s14、计算目标区域内各种土地利用类型对应的单位面积碳储量，包括：
[0149]
步骤s141、在目标区域内设置采样点，获取各个采样点处的各种土地利用类型的单位面积碳储量数据；
[0150]
步骤s142、根据各采样点处的各种土地利用类型的单位面积碳储量数据，通过统计学法获得栅格尺度各种土地利用类型的单位面积碳储量数据，作为第四项指标参数，用于确保生态安全。
[0151]
具体的，步骤s141中，搜集现有研究中关于目标区域上不同采样点处的各种土地利用类型对应的土壤单位面积碳储量、地下部分单位面积碳储量和地上部分单位面积碳储量数据等采样结果。其中，已有研究中采样点的设定是根据目标土地的利用类型、地形、土壤等综合选取。
[0152]
一个优选的实施方式，步骤s142中通过空间插值方法选择广义线性模型，借助r语言编程软件，利用每个采样点的地理位置、采样结果、土地利用类型地形、土壤等参数构建与单位面积碳储量的关系，模拟未知点的单位面积碳储量数据，从而实现对栅格尺度不同
土地利用类型对应土壤单位面积碳储量、地下部分单位面积碳储量和地上部分单位面积碳储量数据的获取，并将其相加得到总的单位面积碳储量数据，作为第四项指标参数，以保证生态安全。
[0153]
其它实施方式中，还可以选择反距离权重插值、克里金插值、样条插值等方法模拟未知点的单位面积碳储量数据。
[0154]
步骤s15、计算目标区域内土地利用类型转变后对应的单位面积物种丰富度，包括：
[0155]
步骤s151、获取基于当前年份土地利用类型分布情况下，栅格尺度的单位面积总物种丰富度数据，并作为单位面积物种丰富度初始数值；
[0156]
步骤s152、计算各种土地利用类型转变后对应的单位面积物种丰富度变化百分比；
[0157]
步骤s153、基于步骤s151和步骤s152的结果计算当前单位面积物种丰富度初始值转变为其他土地利用类型后的单位面积物种丰富度。
[0158]
具体的，步骤s151中，基于世界自然保护联盟(iucn)提供的鸟类、两栖和哺乳类动物的栅格尺度单位面积物种丰富度数据相加，得到栅格尺度单位面积总物种丰富度数据，空间分辨率为10km。将得到的数据作为基于当前土地利用类型分布下，单位面积物种丰富度初始数值。
[0159]
步骤s152中，通过文献资料，计算不同土地利用类型之间转换带来的单位面积物种丰富度变化百分比，以耕地转为林地为例，计算单位面积物种丰富度增加或者减少的百分比。
[0160]
步骤s153中，基于步骤s151计算得到的当前单位面积物种丰富度初始值以及步骤s252得到的不同土地利用类型转换带来的单位面积物种丰富度变化百分比，两个数值相乘即可得到当前土地利用类型转变为其他土地利用类型后的单位面积物种丰富度。该数值作为第五项输入参数，用于保证生态安全。
[0161]
步骤s2、设置优化方法的约束条件，约束条件包括：
[0162]
一、优化后目标区域内每种类型作物的总热量不减少；
[0163]
二、优化后目标区域内总物种丰富度不下降；
[0164]
三、优化后目标区域内总碳储量不减少。
[0165]
通过设置上述三个约束条件，确保优化后的土地布局能够兼顾粮食安全和生态安全。
[0166]
在一个优选实施例中，还包括：
[0167]
四、优化后目标区域内每种类型作物的种植面积不超过该类型作物的种植适宜性面积；
[0168]
五、优化后目标区域内的耕地面积不低于目标耕地面积；
[0169]
六、优化后目标区域内的建设用地面积占比不低于目标建设用地面积占比；
[0170]
七、优化后目标区域内的自然保护区的总物种丰富度不减少；
[0171]
八、优化后目标区域内的水域面积不变；
[0172]
九、优化前未利用地面积占比超过80％或者海拔超过4500m的目标区域，优化后该目标区域内未利用地的面积不变。
[0173]
进一步的，约束条件一中设定总热量不减少是为了保证粮食安全，可以表征为下式：
[0174][0175]
其中，i表示不同栅格单元，i＝1
…
n；
[0176]
j对应于不同土地利用类型，j＝1
…
6；
[0177]
1为耕地、2为林地、3为草地、4为水域、5为建设用地、6为未利用地；
[0178]
l表示作物类型；
[0179]
cal
i，l
表示栅格单元i上第j类土地类型上每100g作物l可食用部分所含卡路里；
[0180]yi，l
表示栅格单元i上第j类土地类型上作物l的单产；
[0181]mi，l
表示栅格单元i上作物l的种植面积占耕地面积的比例；
[0182]
x
i，j
表示优化前栅格单元i上第j类土地类型的用地面积；
[0183]
约束条件二用于保证生态安全，可以表征为下式：
[0184][0185]
其中，o
i，j
表示优化后栅格单元i上第j类土地类型的用地面积；
[0186]bi，j
表示栅格单元i上第j类土地类型对应的单位面积物种丰富度；
[0187]
x
i，j
表示优化前栅格单元i上第j类土地利用类型的用地面积；
[0188]
上式中，单位面积物种丰富度与用地面积的乘积结果表示该目标区域内总物种丰富度。
[0189]
约束条件三用于保证生态安全，可以表征为下式：
[0190][0191]
其中，
[0192]ci，j
表示栅格单元i上第j类土地利用类型对应的单位面积碳储量；
[0193]
上式中，单位面积碳储量与用地面积的乘积结果表示该目标区域内总碳储量。
[0194]
约束条件四、用于约束一些不适合作物种植的地方在优化后有种植，可以表征为下式：
[0195]
0≤m
i，l
×oi，j＝1
≤suitarea
i，l
ꢀꢀ
(7)
[0196]
其中，suitarea
i，l
表示栅格单元i上作物l的种植适宜性面积；
[0197]mi，l
表示栅格单元i上作物l的种植面积占耕地面积的比例。
[0198]
约束条件五、用于防止优化后耕地大幅度减少，不满足总的耕地要求，可以表征为下式：
[0199][0200]
其中，crop表示不同目标区域的耕地目标面积；
[0201]
p表示目标区域内栅格单元i所属的范围；
[0202]
n表示栅格单元的数量。
[0203]
进一步的，通过获取目标区域内当前耕地面积占全国当前总耕地面积的比例，再根据上述比例和全国目标总耕地面积获取目标区域的目标耕地面积。
[0204]
约束条件六、用于防止因水域减少导致缺水情况的发生，可以表征为下式：
[0205]oi，j＝4
＝x
i，j＝4
ꢀꢀ
(9)
[0206]
约束条件七、用于保证建设用地能够满足未来人口增加带来的居住需求，可以表征为下式：
[0207][0208]
其中，built表示不同目标区域的建设用地目标占比。
[0209]
进一步的，通过获取目标区域在不同时段建设用地占土地利用总面积的比例，将不同时段的建设用地占比的平均值作为目标区域的目标建设用地面积占比。
[0210]
一个优选的实施方式，可以结合已有研究模拟的未来不同社会经济发展时段下的土地利用预测数据，借助arcgis软件的“分区统计”工具，计算不同时段下各目标区域的建设用地占该目标区域总面积的比例，取不同时段计算得到比例的平均值作为最终该目标区域内目标建设用地面积占比；与目标区域总面积相乘，得到目标区域的目标建设用地面积。
[0211]
约束条件八的设立是基于原则上未利用地是不能增加的，即未利用地面积占比在80％以上的或者海拔在4500米以上的地区为荒漠或者高寒地区，其土地类型的属性无法转化，可以表征为下式：
[0212]oi，j＝6
＝x
i，j＝6
ꢀꢀ
(10)
[0213]
约束条件九的设立是基于自然保护区设立的目的，即自然保护区本身要保证生态安全，因此生态保护区内的总物种丰富度不能减少，可以表征为下式：
[0214][0215]
其中，z表示自然保护区；
[0216]bi,j
表示栅格单元i上第j类土地类型对应的单位面积物种丰富度。
[0217]
因此，本发明的土地布局优化方法通过约束作物的总热量不减少、物种丰富度不下降以及碳储量不减少以保证粮食安全目标和生态安全目标，同时约束未来建设用地分布、未来建设用地需求和未来耕地需求等，以实现对耕地、建设用地和生态用地等多地类进行综合调整。
[0218]
步骤s3、基于目标函数和约束条件，对目标土地布局进行优化。
[0219]
本发明的土地布局优化方法，基于步骤s1得到的数据，借助地理信息系统软件和通用代数建模系统建立栅格尺度的土地利用布局优化模型，通过对模型设置约束条件、目标函数，优化模拟保证粮食安全和生态安全的各种土地利用布局。
[0220]
因此，在一个优选的实施方式中，还包括利用上述步骤中获得的指标参数，使用软件系统进行建模，实现土地布局优化的步骤，如图3所示：
[0221]
步骤s31、通过arcgis的“创建渔网”工具，创建研究区范围内10km
×
10km距离的渔网点，并提取各渔网点步骤s21-步骤s25得到的各类数据，并保存为excel格式的数据，用于后续优化模拟。
[0222]
步骤s32、在通用代数建模系统(gams)软件中新建工程文件，读取优化算法，选择conopt4求解器，并依据数据自身大小，设定程序运行的内存限制，运行优化算法。
[0223]
步骤s33、将步骤s32优化得到的结果，利用arcgis软件属性表连接工具，连接到步骤s31创建的渔网点上，然后通过“点转栅格”工具，即可得到栅格结果，实现空间化。
[0224]
本实施方式中，优化算法包含：步骤s2中的目标函数和步骤s1中的约束条件两部分。
[0225]
目标函数的设置为保证优化后同时实现作物总热量最大、总物种丰富度最大和总碳储量最大，可以表征为：
[0226][0227]
其中，i表示不同的栅格单元，n表示目标区域栅格单元的总数，i＝1
…
n；
[0228]
j对应于不同土地利用类型，j＝1
…
6；
[0229]
l表示作物类型；
[0230]
w1表示作物总热量的权重；
[0231]
w2表示总碳储量的权重；
[0232]
w3表示总物种丰富度的权重，w1+w2+w3＝1；
[0233]yi，l
表示栅格单元i上作物l的单产；
[0234]
cal
i，l
表示栅格单元i上每100g作物l可食用部分所含卡路里；
[0235]mi，l
表示栅格单元i上作物l的种植面积占耕地面积的比例；
[0236]oi，j
表示优化后栅格单元i上第j类土地类型的用地面积；
[0237]
x
i，j
表示优化前栅格单元i上第j类土地类型的用地面积；
[0238]ci，j
表示栅格单元i上第j类土地类型对应的单位面积碳储量；
[0239]bi，j
表示栅格单元i上第j类土地类型对应的单位面积物种丰富度；
[0240]
本发明的土地布局优化方法通过约束作物的总热量不减少、总物种丰富度不下降以及总碳储量不减少以保证粮食安全目标和生态安全目标，同时约束未来建设用地分布、未来建设用地需求和未来耕地需求等，以实现对耕地、建设用地和生态用地等多地类进行综合调整。一方面通过约束条件的限制，保证了三者的基本安全；另一方面通过目标函数的设定，保证三者同时达到最优水平。有助于国家或地方政府更加方便的展开有针对性的、具体的国土空间规划措施，从栅格单元进行管理、调整，从而有效的协调耕地、建设用地和生态用地之间的竞争，以实现粮食安全和生态安全。
[0241]
实施例二：
[0242]
中国作为全球生物多样性最丰富的12个国家之一，拥有森林、草原、荒漠、湿地等重要的陆地生态系统，并拥有34984种高等植物，6445种脊椎动物。但目前我国生物多样性下降的总体趋势尚未得到有效遏制，部分生态系统功能不断退化，物种濒危程度加剧。此外，我国作为农业大国，农田扩张等通过导致栖息地丧失威胁生物多样性，并通过生物量和土壤碳的损失影响碳储存。我国城镇化率的不断增加，也进一步影响生物多样性等。如何在保障我国未来粮食需求和建设用地需求基础上，降低对生物多样性和碳储量的影响成为研究重点。
[0243]
因此，本发明第二实施例以中国为目标区域，对2015年中国的土地利用布局进行
优化调整，进一步阐述本方法的可行性和有效性。
[0244]
下面对本实施例中的研究数据来源进行说明：
[0245]
2015年中国土地利用百分比栅格数据来源于中国科学院资源环境科学与数据中心；
[0246]
2015年全国各县区尺度作物播种面积和单产数据来源于各类农业统计年鉴；
[0247]
2010年栅格尺度作物的播种面积和单产来源于spam2010数据集；
[0248]
陆地生态系统的单位面积碳储量数据来源于国家生态科学数据中心；
[0249]
物种丰富度数据来源于世界自然保护联盟(iucn)；
[0250]
各种土地利用类型下的物种丰富度百分比数据来源于projects数据库；
[0251]
不同时段的建设用地的模拟数据来源于geosos网站；
[0252]
作物熟制数据来源于已有研究数据集。
[0253]
步骤s1、获取目标函数的输入参数，包括下述步骤：
[0254]
步骤s11、基于土地利用百分比栅格数据，计算目标区域上每个栅格单元内不同土地利用类型的用地面积；
[0255]
具体的，将中国土地划分成若干个空间分辨率为10km的栅格单元，利用2015年中国土地利用百分比栅格数据，即每个栅格单元内不同土地利用类型的面积占该栅格单元总面积的比例，可计算得到每个栅格单元内各种土地利用类型的面积。
[0256]
本实施例中，土地利用类型包括：耕地、建设用地、林地、草地、水域和未利用地共六大类，每个栅格单元内均可包含上述六类土地利用类型。由于空间分辨率为10km，栅格单元总面积为100km2，即10000ha。
[0257]
步骤s12、计算目标区域内各栅格单元内各种类型作物的作物播种面积和作物单产。
[0258]
具体的，搜集2015年中国各县区小麦、玉米和水稻的播种面积和单产数据作为第一作物产量数据集，以及公开发表的spam2010数据集spam(spatial production allocation model)作为第二作物产量数据集。
[0259]
结合spam2010数据集中的栅格尺度播种面积和单产，将2015年的县区级尺度的作物数据降尺度为栅格尺度的作物数据，得到2015年栅格尺度的作物播种面积和单产。栅格单元分辨率为10km
×
10km。
[0260]
进一步的，作物播种面积的分配算法为：
[0261][0262]
其中，r为栅格单元的编号；c为县区；
[0263]
f为作物类型；y为2015年；
[0264]
ha
r，c，f，y
表示2015年中，某县区c的栅格单元r上作物f的播种面积；
[0265]
ha
r，c，f，2010
表示第二作物产量数据集中，某县区c的栅格单元r上作物f的播种面积；
[0266]
staha
c，f，y
表示第一作物产量数据集中，2015年中，某县区c栅格单元r上作物f的播种面积。
[0267]
作物单产的分配算法为：
[0268][0269][0270]
其中，y
r，c，f，2010
表示第二作物产量数据集中某县区c上栅格单元r上作物f的单产；
[0271]
yavg
c，f，2010
表示由第二作物产量数据集计算得到的某县区c上作物f的平均单产；
[0272]
stay
c，f，y
表示第一作物产量数据集中2015年，某县区c上作物f的平均单产；
[0273]yr,c,f,y
表示2015年，某县区c上栅格单元r上作物f的单产。
[0274]
步骤s13、计算每个栅格单元内作物的种植面积与耕地面积的比例，
[0275]
具体的，搜集栅格尺度作物熟制数据，结合步骤s12计算得到的栅格尺度作物播种面积，计算得到栅格尺度作物的种植面积，即作物种植实际占用耕地的面积。播种面积是考虑作物熟制后的作物总种植面积，通过与作物熟制相除，即可得到作物种植面积。再通过与步骤s11计算得到的每个栅格单元内的耕地面积相除，即可得到作物种植面积占耕地的比例。
[0276]
步骤s14、计算目标区域栅格尺度的不同土地利用类型的单位面积碳储量；
[0277]
具体的，首先，通过查阅相关文献，搜集全国范围内不同采样点的各类土地利用类型对应的土壤单位面积碳储量、地下部分单位面积碳储量和地上部分单位面积碳储量数据，并通过空间插值的方法得到栅格尺度不同土地利用类型对应土壤单位面积碳储量、地下部分单位面积碳储量和地上部分单位面积碳储量数据，并将其相加，得到总的单位面积碳储量数据。空间插值方法选择广义线性模型，借助r语言编程软件实现。
[0278]
步骤s15、计算目标区域内各栅格单元内各种土地利用类型的单位面积物种丰富度；
[0279]
具体的，包括步骤：
[0280]
步骤s151、基于世界自然保护联盟(iucn)提供的鸟类、两栖和哺乳类动物的栅格尺度单位面积物种丰富度数据相加，得到栅格尺度单位面积总物种丰富度数据，空间分辨率为10km。得到的数据表示当前的分布情况和初始数值；
[0281]
步骤s152、通过文献资料，计算不同土地利用类型之间转换带来的单位面积物种丰富度变化百分比，如耕地转为林地，单位面积物种丰富度增加或者减少的百分比；
[0282]
步骤s153、基于步骤s151计算得到的2015年单位面积物种丰富度初始值以及步骤s252得到的不同土地利用类型转换带来的单位面积物种丰富度变化百分比，即可计算2015年度物种丰富度初始值转变为其他土地利用类型后的单位面积物种丰富度。
[0283]
步骤s16、计算各省份的目标耕地面积和目标建设用地面积占比，包括：
[0284]
步骤s161、计算当前各省耕地面积占全国总耕地面积的比例，将中国18.25亿亩耕地保护红线目标按比例计算分配，得到18.25亿亩目标下各省目标耕地面积；
[0285]
步骤s162、结合已有研究模拟的未来不同社会经济发展情景下的土地利用预测数据，借助arcgis软件的“分区统计”工具，计算不同社会经济发展时段下各省建设用地占各省总面积的比例，取不同情景计算得到比例的平均值作为最终的目标建设用地面积占比。
[0286]
步骤s2、设置下述约束条件：
[0287]
一、优化后全国范围内每种类型作物的总热量不减少；
[0288]
二、优化后全国范围内总物种丰富度不下降；
[0289]
三、优化后全国范围内总碳储量不减少；
[0290]
四、优化后目标区域内各栅格单元内每种类型作物的种植面积不超过该栅格单元内该类型作物的种植适宜性面积；
[0291]
五、优化后各省份的耕地面积不低于各省份的目标耕地面积；
[0292]
六、优化后各省份的建设用地面积占比不低于各省份的目标建设用地面积占比；
[0293]
七、优化后自然保护区内每个栅格内的总物种丰富度不减少；
[0294]
八、优化后每个栅格单元内的水域面积不变；
[0295]
九、优化前未利用地面积占比超过80％或者海拔超过4500m的目标区域，优化后该目标区域内未利用地的面积不变。
[0296]
步骤s3、基于目标函数和约束条件，对目标土地利用布局进行优化。
[0297]
基于步骤s1得到的数据，借助地理信息系统软件和通用代数建模系统建立栅格尺度的土地利用布局优化模型，通过对模型设置目标函数、约束条件，优化模拟保证粮食安全和生态安全的土地利用布局，包括：
[0298]
步骤s31、通过arcgis的“创建渔网”工具，创建研究区范围内10km
×
10km距离的渔网点，并提取各渔网点步骤s21-步骤s25得到的各类数据，并保存为excel格式的数据，用于后续优化模拟。
[0299]
步骤s32、在通用代数建模系统(gams)软件中新建工程文件，读取优化算法，选择conopt4求解器，并依据数据自身大小，设定程序运行的内存限制，运行优化算法。
[0300]
步骤s33、将步骤s32优化得到的结果，利用arcgis软件属性表连接工具，连接到步骤s31创建的渔网点上，然后通过“点转栅格”工具，即可得到栅格结果，实现空间化。得到如下研究结果：
[0301]
(1)优化调整后耕地、建设用地和生态用地变化分析：
[0302]
应用本优化算法，通过对耕地、建设用地以及生态用地的综合调整，使得粮食安全和生态安全同时实现。优化结果显示，在空间上，耕地仍然主要分布在东部，尤其是东北平原、黄淮海平原以及长江中下游地区，这些地区耕地面积也以增加为主，主要是由于这些地区作物单产水平高，为满足粮食安全目标，种植面积有所增加。生态用地整体仍然以西部和南部分布最为广泛，增加地区也主要集中在南部地区。主要是由于南部地区林草地的物种丰富度和碳储量相对较高，为满足生态安全目标，南部地区呈现增加趋势。而建设用地则呈现零星分布，主要是满足未来建设用地需求。从具体数值上来看，相对于原布局，优化调整后全国耕地增加了9.49％，建设用地增加了18.23％，生态用地减少了2.92％。
[0303]
(2)优化调整后作物总热量、总物种丰富度和总碳储量变化分析
[0304]
优化调整后，不仅能够保证作物总热量、总物种丰富度和总碳储量不减少，还能最大程度的使得三者最大化，实现粮食、生态的绝对安全。结果统计显示，作物总热量增加了32.80％，其中小麦总热量增加18.26％，玉米总热量增加42.84％，水稻总热量增加28.76％。而相应的总物种丰富度增加5.72％，总碳储量增加15.64％。
[0305]
应当理解的是，本发明的上述具体实施方式仅仅用于示例性说明或解释本发明的原理，而不构成对本发明的限制。因此，在不偏离本发明的精神和范围的情况下所做的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。此外，本发明所附权利要求旨
在涵盖落入所附权利要求范围和边界、或者这种范围和边界的等同形式内的全部变化和修改例。

技术特征：
1.一种土地布局优化方法，其特征在于，包括：设置优化方法的目标函数，所述目标函数的优化目标至少包括：实现作物总热量最大、总物种丰富度最大和总碳储量最大；设置优化方法的约束条件，所述约束条件至少包括：优化后目标区域内每种类型作物的总热量不减少、优化后目标区域内总物种丰富度不下降，以及优化后目标区域内总碳储量不减少；基于所述约束条件和目标函数，对目标土地布局进行优化。2.根据权利要求1所述的土地布局优化方法，其特征在于，所述约束条件还包括：优化后目标区域内每种类型作物的种植面积不超过该类型作物的种植适宜性面积；优化后目标区域内的耕地面积不低于目标耕地面积；优化后目标区域内的建设用地面积占比不低于目标建设用地面积占比；优化后目标区域内的自然保护区的总物种丰富度不减少；优化后目标区域内的水域面积不变；优化前未利用地面积占比超过80％或者海拔超过4500m的目标区域，优化后该目标区域内未利用地的面积不变。3.根据权利要求1或2所述的土地布局优化方法，其特征在于，所述目标函数的输入参数包括：目标区域内各种土地利用类型的用地面积；目标区域内各种类型作物的作物播种面积和作物单产；目标区域内作物种植面积占用耕地的比例；目标区域内各种土地利用类型对应的单位面积碳储量；目标区域内土地利用类型转变后对应的单位面积物种丰富度。4.根据权利要求3所述的土地布局优化方法，其特征在于，所述目标区域内各种土地利用类型的用地面积，通过下述步骤获取：获取所述目标区域的土地利用百分比栅格数据；利用所述土地利用百分比栅格数据，计算目标区域上每个栅格单元内各种土地利用类型的用地面积。5.根据权利要求3所述的土地布局优化方法，其特征在于，所述目标区域内各种类型作物的作物播种面积和作物单产，通过下述步骤获取：获取目标区域的各种类型作物的作物产量数据集；计算目标区域上某年份每个栅格单元内各种类型作物的作物播种面积；计算目标区域上某年份每个栅格单元内各种类型作物的作物单产；所述作物产量数据集包括：第一作物产量数据集和第二作物产量数据集；所述第一作物产量数据集为具有各种作物类型的作物播种面积和作物单产的统计数据集；所述第二作物产量数据集为具有各种作物类型的栅格尺度的作物播种面积和作物单产的数据集。6.根据权利要求5所述的土地布局优化方法，其特征在于，所述作物播种面积的分配算法为：
其中，r为栅格单元的编号；c为目标区域；f为作物类型；y为模拟年份；ha
r,c,f,y
表示模拟年份y中，目标区域c上的栅格单元r上作物f的播种面积；ha
r,c,f,2010
表示第二作物产量数据集中，目标区域c上栅格单元r上作物f的播种面积；staha
c,f,y
表示第一作物产量数据集中，模拟年份y中，目标区域c上栅格单元r上作物f的播种面积；所述作物单产的分配算法为：所述作物单产的分配算法为：其中，y
r,c,f,2010
表示第二作物产量数据集中目标区域c上栅格单元r上作物f的单产；yavg
c,f,2010
表示由第二作物产量数据集计算得到的目标区域c上作物f的平均单产；stay
c,f,y
表示第一作物产量数据集中模拟年份y中，目标区域c上作物f的平均单产；y
r,c,f,y
表示模拟年份y中，目标区域c的栅格单元r上作物f的单产。7.根据权利要求3所述的土地布局优化方法，其特征在于，所述各种土地利用类型对应的单位面积碳储量，通过下述步骤获取：在所述目标区域内设置采样点，获取各个采样点处的各种土地利用类型的单位面积碳储量数据；根据所述各采样点处的各种土地利用类型的单位面积碳储量数据，利用统计学方法获得栅格尺度各种土地利用类型的单位面积碳储量数据。8.根据权利要求3所述的土地布局优化方法，其特征在于，所述目标区域内土地利用类型转变后对应的单位面积物种丰富度，通过下述步骤获取：获取基于当前年份土地分布情况下栅格尺度的单位面积总物种丰富度数据，并作为单位面积物种丰富度初始数值；计算各种土地利用类型之间转换带来的单位面积物种丰富度变化百分比；基于所述单位面积物种丰富度初始数值和所述单位面积物种丰富度变化百分比，计算得到变化后的目标区域内土地利用类型转变后对应的单位面积物种丰富度。9.根据权利要求2所述的土地布局优化方法，其特征在于，所述目标耕地面积和所述目标建设用地面积，通过下述步骤获取：获取目标区域内当前耕地面积占全国当前总耕地面积的比例；根据上述比例和全国目标总耕地面积获取目标区域的目标耕地面积；获取目标区域在不同时段的建设用地占土地利用总面积的比例；将不同时段的建设用地占比的平均值作为目标区域的目标建设用地面积占比。10.根据权利要求1-9任一项所述的土地布局优化方法，其特征在于，所述目标函数表征为：
其中，i表示不同的栅格单元，n表示目标区域栅格单元的总数，i＝1
…
n；j表示不同土地利用类型，j＝1
…
6；l表示作物类型；w1表示作物总热量的权重；w2表示总碳储量的权重；w3表示总物种丰富度的权重，w1+w2+w3＝1；y
i,l
表示栅格单元i上作物l的单产；cal
i,l
表示栅格单元i上每100g作物l可食用部分所含卡路里；m
i,l
表示栅格单元i上作物l的种植面积占耕地面积的比例；o
i,j
表示优化后栅格单元i上第j类土地类型的用地面积；x
i,j
表示优化前栅格单元i上第j类土地类型的用地面积；c
i,j
表示栅格单元i上第j类土地类型对应的单位面积碳储量；b
i,j
表示栅格单元i上第j类土地类型对应的单位面积物种丰富度。

技术总结
本发明公开了一种土地布局优化方法，其包括步骤：设置优化方法的目标函数；设置优化方法的约束条件；基于所述目标函数和约束条件，对目标土地布局进行优化。本发明的土地布局优化方法在栅格尺度上通过优化不同土地利用类型分布，通过约束作物的总热量不减少、物种丰富度不下降以及碳储量不减少以保证粮食安全目标和生态安全目标，同时约束未来建设用地分布、未来建设用地需求和未来耕地需求等，以实现对耕地、建设用地和生态用地等多地类进行综合调整。合调整。合调整。


技术研发人员：史文娇 王鸣雷 岳天祥
受保护的技术使用者：中国科学院地理科学与资源研究所
技术研发日：2023.06.05
技术公布日：2023/9/9